CACCS: Secuencia de taller de RStudio – Parte 3
Rashid C.J. Marcano Rivera
10 de oct.ᵉ de 2025
Inferencia estadística
Este taller está basado en elementos y ejemplos del libro de Rafael Irizarry, aquí, así como pedazos del curso de la Universidad de la República, en Uruguay, disponible en esta página de RPubs de RStudio y este repaso sobre modelos longitudinales por Alessio Crippa también en RPubs de RStudio. Recomiendo complementar el análisis con el libro Data Analysis Using Regression and Multilevel/Hierarchical Models de Andrew Gelman y Jennifer Hill.
Si aún no has instalado R, está aquí. Acto seguido, baja RStudio. Puedes también ir a la nube en Posit Cloud.
1 Recapitulando
La vez anterior, tomamos un recorrido a través de distintos tipos de visualizaciones. En efecto, una buena visualización puede demostrar bastante sobre elementos en nuestros estudios. Por ejemplo
library(wooldridge)
data(wage1)
head(wage1)## wage educ exper tenure nonwhite female married numdep smsa northcen south
## 1 3.10 11 2 0 0 1 0 2 1 0 0
## 2 3.24 12 22 2 0 1 1 3 1 0 0
## 3 3.00 11 2 0 0 0 0 2 0 0 0
## 4 6.00 8 44 28 0 0 1 0 1 0 0
## 5 5.30 12 7 2 0 0 1 1 0 0 0
## 6 8.75 16 9 8 0 0 1 0 1 0 0
## west construc ndurman trcommpu trade services profserv profocc clerocc
## 1 1 0 0 0 0 0 0 0 0
## 2 1 0 0 0 0 1 0 0 0
## 3 1 0 0 0 1 0 0 0 0
## 4 1 0 0 0 0 0 0 0 1
## 5 1 0 0 0 0 0 0 0 0
## 6 1 0 0 0 0 0 1 1 0
## servocc lwage expersq tenursq
## 1 0 1.131402 4 0
## 2 1 1.175573 484 4
## 3 0 1.098612 4 0
## 4 0 1.791759 1936 784
## 5 0 1.667707 49 4
## 6 0 2.169054 81 64¿Qué aprendemos de ver estos datos así? ¿Podemos rápidamente determinar a si años de educación se traducen a mayores ingresos? ¿Podemos determinar si afecta, en algo, la relación marital? Para muchos humanos, es difícil extraer información con meramente mirar a números sin contexto adicional. Pero podríamos ver algo en este gráfico
Vivimos en una era de creciente disponibilidad de conjuntos de datos informativos y de herramientas de software, con lo cual el uso de visualizaciones ha aumentado en diversos espacios: académicos, gubernamentales, organizaciones sociales, prensa, e industrias varias. Sin embargo, R, un programa estadístico diseñado para el manejo y análisis de distintos formatos de datos, tenemos un sinnúmero de opciones para trabajar distintos tipos de complejidades en datos y análisis.
En esta secuencia de cuatro semanas, continuamos de lo aprendido en el pasado; aplicaremos funciones para cargar datos, manipularlos, y retomamos problemas que visualizamos la semana pasada como el señalado, para analizarlos con análisis estadísticos y visualizaciones relativas a la inferencia estadística. De una vez estaremos entrando en algunos de los diagnósticos que podremos desplegar para evaluar qué tan adecuado es un modelo para trabajar ciertos datos. Aprenderemos hoy:
Varias operaciones estadísticas,
Estadística inferencial, en sus varias versiones para modelos simples, lineales, jerárquicos y longitudinales.
Gráficos y pruebas que ayuden a entender y diagnosticar los modelos estadísticos que usaremos.
2 Modelos estadísticos
Queríamos entender la vez anterior la relación de ingreso con otras variables. Para cargar los datos escribiremos
library(wooldridge)
base <- wage1
#View(base)
names(base)## [1] "wage" "educ" "exper" "tenure" "nonwhite" "female"
## [7] "married" "numdep" "smsa" "northcen" "south" "west"
## [13] "construc" "ndurman" "trcommpu" "trade" "services" "profserv"
## [19] "profocc" "clerocc" "servocc" "lwage" "expersq" "tenursq"#?wage1Las variables que utilizaremos son las siguientes:
wage: salario promedio por hora.
educ: años de educación.
exper: años de experiencia potencial.
tenure: años con el empleador actual (antigüedad).
nonwhite: es igual a 1 si la persona no es blanca, 0 si no.
female: es igual a 1 si la persona es mujer, 0 si no
married: es igual a 1 si la persona es casada, 0 si no.
En primer lugar queremos cambiar el nombre de la variable que está en la posición 4:
names(base)[4] <- "antigüedad"
names(base)[24] <- "antigüedadcuad"Por ahora lo que nos interesaba es un subconjunto de variables, todas de la 1 a la 7 (la 4.ª ha quedado como antigüedad), la 22 (log. de salario), la 23 (experiencia al cuadrado) y la 24 (la antigüedad al cuadrado).
base1 <- base[,c(1: 7, 22:24)] Veremos ahora las primeras filas
head(base1, n = 10)## wage educ exper antigüedad nonwhite female married lwage expersq
## 1 3.10 11 2 0 0 1 0 1.131402 4
## 2 3.24 12 22 2 0 1 1 1.175573 484
## 3 3.00 11 2 0 0 0 0 1.098612 4
## 4 6.00 8 44 28 0 0 1 1.791759 1936
## 5 5.30 12 7 2 0 0 1 1.667707 49
## 6 8.75 16 9 8 0 0 1 2.169054 81
## 7 11.25 18 15 7 0 0 0 2.420368 225
## 8 5.00 12 5 3 0 1 0 1.609438 25
## 9 3.60 12 26 4 0 1 0 1.280934 676
## 10 18.18 17 22 21 0 0 1 2.900322 484
## antigüedadcuad
## 1 0
## 2 4
## 3 0
## 4 784
## 5 4
## 6 64
## 7 49
## 8 9
## 9 16
## 10 441También podríamos llamarlos con el nombre de variable
datos1 <- base[,c("wage","educ","exper", "antigüedad" )]
head(datos1, n = 5)## wage educ exper antigüedad
## 1 3.10 11 2 0
## 2 3.24 12 22 2
## 3 3.00 11 2 0
## 4 6.00 8 44 28
## 5 5.30 12 7 2o como cubrimos al tocar tidyverse, la función
select()
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.1 ✔ stringr 1.5.2
## ✔ ggplot2 4.0.0 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorscolnames(base)## [1] "wage" "educ" "exper" "antigüedad"
## [5] "nonwhite" "female" "married" "numdep"
## [9] "smsa" "northcen" "south" "west"
## [13] "construc" "ndurman" "trcommpu" "trade"
## [17] "services" "profserv" "profocc" "clerocc"
## [21] "servocc" "lwage" "expersq" "antigüedadcuad"datos2 <- base |>
dplyr::select(wage, educ, exper, antigüedad) #¿qué pasa si no tengo la especificación?
head(datos2, n=7)## wage educ exper antigüedad
## 1 3.10 11 2 0
## 2 3.24 12 22 2
## 3 3.00 11 2 0
## 4 6.00 8 44 28
## 5 5.30 12 7 2
## 6 8.75 16 9 8
## 7 11.25 18 15 7Si solamente quisiéramos los datos de los casados podríamos usar la función de filtrado
datos3 <- base |>
dplyr::select("wage", "educ", "exper", "antigüedad", "married") |>
filter(married == 1)3 Correlaciones
Para investigar si hay correlación entre alguna de las variables se puede realizar un gráfico en el que se presenta la dispersión para cada par de variables.
plot(datos2)
Y también calcular la matriz de correlaciones de las variables que figuran en datos2.
cor(datos2)## wage educ exper antigüedad
## wage 1.0000000 0.40590333 0.1129034 0.34688957
## educ 0.4059033 1.00000000 -0.2995418 -0.05617257
## exper 0.1129034 -0.29954184 1.0000000 0.49929145
## antigüedad 0.3468896 -0.05617257 0.4992914 1.00000000Notamos que existe una correlación positiva entre el salario y la educación (0.4059)
Crearemos el diagrama de dispersión entre salario y educación utilizando las funciones de la librería ggplot2.
ggplot(datos2, aes(x = educ, y = wage)) +
geom_point() + theme_light() +
ggtitle("Relación entre salario y educación")
Se acordarán que la vez anterior añadimos complicaciones como un nivel adicional en la capa de color:
wage1 |>
mutate(marital = factor(married, levels = c(0, 1), labels = c("Soltero", "Casado"))) |>
ggplot(aes(educ, wage)) + geom_point(aes(colour = marital), size = 3)+
labs(title="Relación entre educación y salario",
x = "Educación",
y = "Salario",
color = "Estado marital",
caption = "Datos de Wooldridge, usado en págs. como 7, 17, 33-34, 37, 76,...")
Esto lo exploraremos más al entrar en modelos de regresión lineal múltiple. Empezaremos por el modelo de regresión lineal simple.
4 Modelo de regresión lineal simple
Estimamos un modelo de regresión lineal simple, con el método mínimos
cuadrados ordinarios (en adelante MCO, OLS en inglés) que explique los
salarios en función de los años de educación de las personas. En R esto
se hace con la función lm():
mod1 <- lm(wage ~ educ, data=datos2)
summary(mod1) # para imprimir la salida##
## Call:
## lm(formula = wage ~ educ, data = datos2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.3396 -2.1501 -0.9674 1.1921 16.6085
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.90485 0.68497 -1.321 0.187
## educ 0.54136 0.05325 10.167 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.378 on 524 degrees of freedom
## Multiple R-squared: 0.1648, Adjusted R-squared: 0.1632
## F-statistic: 103.4 on 1 and 524 DF, p-value: < 2.2e-16¿Cómo se lee esta salida? La primera línea indica la fórmula que se
utilizó. La segunda es sobre la distribución de residuos las diferencias
entre los valores observados de wage y los valores
predichos por el modelo. La mediana cercana a cero indica que los
residuos están centrados alrededor de cero. El rango de los residuos
sugiere que hay algunos valores atípicos, especialmente en el extremo
máximo (16.6085), lo que podría indicar la presencia de salarios
excepcionalmente altos no explicados completamente por el modelo.
Vamos a los coeficientes: la información se presenta a través de
varias columnas. La columna de Estimate tiene el valor
estimado de coeficientes (\(\hat{\beta}\), mientras que la columna de
Std. Error nos da en promedio lo que varía el estimado en
relación a la variable dependiente. Esto es de utilidad para computar
intervalos de confianza y establecer la métrica con la cual determinar
la hipótesis la existencia de una relación entre una variable y otra. El
puntaje t reporta la distancia estandarizada en distribución t de
nuestro coeficiente, en relación a la posibilidad de que tuviera cero
efecto. Mientras mayor sea el número de puntaje, y se mantuvieran
relativamente mayores en relación al error estándar indica que una
relación existe. La última columna provee el valor p, probabilidad
asociada al puntaje t, que indica la probabilidad de estar viendo un
valor tan extremo como el reportado por azar.
El intercepto de -0.9 representa el salario promedio cuando los años de educación son cero. Sin embargo, en la práctica, es poco común que una persona tenga cero años de educación (o que tenga ingreso negativo por un trabajo), por lo que este valor tiene una interpretación limitada. El valor p es algo elevado (es decir, supera los niveles utilizados normalmente como corte para evitar errores tipo 1), lo que indica que no podríamos afirmar con certeza su diferencia de cero. Por otro lado, el coeficiente en educación indica que cada año en educación (una unidad adicional) se traduciría en un aumento de 0.54136 unidades en salario (si son pesos, pues 54 chavos). El valor p es bajo (o el puntaje t es elevado, distante a 2, y mucho mayor que el error estándar), lo que indica con cierta certeza que el estimado no es cero. Vemos esto acompañado con asteriscos, simbolizando el nivel de significancia, atado al valor seleccionado \(\alpha\), el punto de corte anteriormente mencionado.
El modelo luego continúa reportando otros diagnósticos generales:
Residual standard error(Error estándar de los residuos): más o menos 3.378- Indica la variabilidad promedio de los residuos; en otras palabras, mide la precisión del modelo. Los modelos lineales incluyen un término estocástico, que captura las desviaciones de nuestra relación predicha con las observaciones. En este caso, el valor de 3.378, sugiere que las predicciones individuales del salario pueden variar en promedio ±3.378 unidades del valor real a través del espacio de nuestras observaciones.
Degrees of freedom(Grados de libertad): 524 – Calculado como el número de observaciones menos el número de parámetros estimados (n - k). El objeto utilizado para análisis tenía 526 filas, y estimamos la pendiente y el intercepto.El coeficiente de determinación (\(R^2\)) nos da 0.1648, que indica que el 16.48% de la variabilidad se explica con el modelo: es decir una proporción de la varianza explicada. El R cuadrado ajustado es similar, pero penaliza la inclusión de términos adicionales, algo que el \(R^2\) no hará: al seguir añadiendo variables, \(R^2\) seguirá aumentando, sean buenas o no las adiciones; el \(R^2\) ajustado es preferido para análisis en regresión múltiple.
El estadístico F y el valor p global evalúa la hipótesis nula de que todos los coeficientes sean iguales a cero.
Al objeto creado para almacenar el modelo lineal le podemos hacer
análisis varios. Por ejemplo, podríamos hacerle el análisis de
varianzas. Este se calcula con la función anova():
anova(mod1) # para imprimir el análisis de varianzas## Analysis of Variance Table
##
## Response: wage
## Df Sum Sq Mean Sq F value Pr(>F)
## educ 1 1179.7 1179.73 103.36 < 2.2e-16 ***
## Residuals 524 5980.7 11.41
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Este análisis descompone la variabilidad total de la variable
dependiente (wage en este caso) en componentes atribuibles
al modelo y a los residuos. ¿Para qué haríamos esto? Quizás en este caso
simple no parezca muy revelador (aunque da más detalles, por ejemplo,
sobre el cómputo que lleva al puntaje F, viendo). Si tenemos variables
categóricas con más de dos niveles, podremos aclarar más el impacto de
la variable en general, en lugar de meramente entender la diferencia de
distintos grupos en referencia a un valor base.
El objeto creado a través de la función lm(), de clase
lm, tiene en sí doce componentes:
names(mod1)## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"Se puede acceder a ellos como a las variables dentro de un objeto,
utilizando el operador $ entre el objeto y el elemento. Por
ejemplo, para extraer los coeficientes estimados se escribe lo
siguiente:
mod1$coefficients## (Intercept) educ
## -0.9048516 0.5413593coef(summary(mod1))## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.9048516 0.68496782 -1.321013 1.870735e-01
## educ 0.5413593 0.05324804 10.166746 2.782599e-22Notemos que la función de coef, que extrae coeficientes
de objetos creados por funciones de modelaje estadístico, nos retorna
una matriz con cuatro columnas y la cantidad de filas adecuadas para los
parámetros estimados, en este caso dos.
Una variación podría ser:
coefficients(mod1)## (Intercept) educ
## -0.9048516 0.5413593summary(mod1)$coefficients #noten, es igual en resultado a la función de coef(summary(modelo))## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.9048516 0.68496782 -1.321013 1.870735e-01
## educ 0.5413593 0.05324804 10.166746 2.782599e-22Para consultar la estimación de un coeficiente de regresión se
utilizan corchetes [ ] y se indica la ubicación del mismo
dentro de la salida del summary(). Por ejemplo, al alfa
gorro o sombrerito (es decir, \(\hat{\alpha}\), alfa estimado) se podría
sacar llamando a la [primera fila, primera columna]:
ahat <- coef(summary(mod1))[1,1]
ahat## [1] -0.9048516Y el beta gorro \(\hat{\beta}\) o beta estimado se podría obtener llamando a la [segunda fila, primera columna]:
bhat <- coef(summary(mod1))[2,1]
bhat## [1] 0.5413593Del modelo podemos sacar también:
coeficientes <- mod1$coefficients #vector de coeficientes estimados
ygorro <- mod1$fitted.values #valores predichos
resid <- mod1$residuals #residuos
#los valores estimados o predichos también se pueden sacar con la función predict():
ygor1<- predict(mod1)
head(ygor1)## 1 2 3 4 5 6
## 5.050100 5.591459 5.050100 3.426022 5.591459 7.756896head(ygorro)## 1 2 3 4 5 6
## 5.050100 5.591459 5.050100 3.426022 5.591459 7.756896all.equal(ygor1, ygorro)## [1] TRUEAgregaremos al conjunto de datos original, las predicciones, en este
caso salario estimado, con la función predict().
datos2$predicciones <- predict(mod1)
head(datos2, 6)## wage educ exper antigüedad predicciones
## 1 3.10 11 2 0 5.050100
## 2 3.24 12 22 2 5.591459
## 3 3.00 11 2 0 5.050100
## 4 6.00 8 44 28 3.426022
## 5 5.30 12 7 2 5.591459
## 6 8.75 16 9 8 7.756896El gráfico de dispersión puede establecerse
ggplot(datos2, aes(x = educ, y = wage)) +
geom_point() +
geom_smooth(method = 'lm', formula = y ~ x, se = FALSE, col = 'dodgerblue1') +
theme_light() +
ggtitle("Relación entre salario y educación")
Si quisiéramos un gráfico de dispersión interactivo, podemos usar plotly. Así, posicionándose encima de cada observación, se ven los valores de (x, y) para cada uno de los individuos. Para construir dicho gráfico se necesita la función ggplotly() del paquete plotly.
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutggplotly(data = datos2, x = ~ educ, y = ~ wage)4.1 Supuestos de modelos lineales
Ahora, a la hora de evaluar si nuestro modelo es bueno, además de los estimados que vimos sobre el ajuste del modelo o su capacidad explicativa de la varianza, tenemos que revisar las asunciones o presunciones que hace un modelo sobre los datos que evalúa. En la evaluación de un modelo lineal, tenemos cinco.
Linealidad: La relación entre las variables dependiente (y) e independiente(s) (x) debe ser lineal en los parámetros. Esto significa que el modelo es lineal en los coeficientes, aunque las variables en sí mismas puedan estar transformadas (por ejemplo, mediante logaritmos).
Independencia de los errores: Los residuos o errores (la diferencia entre los valores observados y los valores predichos por el modelo) deben ser independientes entre sí. Esto significa que no debe haber correlación entre los errores.
Homoscedasticidad: La varianza de los errores debe ser constante a lo largo de todos los valores de las variables independientes. Esto implica que la dispersión de los residuos debe ser más o menos la misma a lo largo del rango de valores de la variable independiente. En caso de heteroscedasticiddad, las estimaciones se tornan ineficientes y los errores estándar incorrectos.
Normalidad de los errores: Los errores deben seguir una distribución normal. Este supuesto es importante para la realización de pruebas de hipótesis y la construcción de intervalos de confianza. Cabe señalar que en este caso, los estimadores de los coeficientes siguen siendo insesgados incluso si los errores no fueran normales. La falta de normalidad afecta principalmente la inferencia estadística.
No multicolinealidad: Las variables independientes no deben estar altamente correlacionadas entre sí. La multicolinealidad puede dificultar la estimación precisa de los coeficientes de regresión.
De violarse estas presunciones, el modelo estará sesgado y sus resultados no serán del todo fiables.
Podremos revisar algunos de estos a través de varios métodos: Podemos por ejemplo calcular los residuos del modelo simple y los agregamos al conjunto de datos (datos2) de la siguiente forma:
datos2$residuos <- datos2$wage - datos2$predicciones
head(datos2, 5)## wage educ exper antigüedad predicciones residuos
## 1 3.10 11 2 0 5.050100 -1.9501003
## 2 3.24 12 22 2 5.591459 -2.3514594
## 3 3.00 11 2 0 5.050100 -2.0501002
## 4 6.00 8 44 28 3.426022 2.5739776
## 5 5.30 12 7 2 5.591459 -0.2914593Esto es lo mismo que R hace internamente, y se puede llamar con la
función residuals:
datos2$residuosmod <- residuals(mod1)
head(datos2, 5)## wage educ exper antigüedad predicciones residuos residuosmod
## 1 3.10 11 2 0 5.050100 -1.9501003 -1.9501003
## 2 3.24 12 22 2 5.591459 -2.3514594 -2.3514594
## 3 3.00 11 2 0 5.050100 -2.0501002 -2.0501002
## 4 6.00 8 44 28 3.426022 2.5739776 2.5739776
## 5 5.30 12 7 2 5.591459 -0.2914593 -0.2914593Comparamos los residuos calculados manualmente con los que nos dio el
modelo. Otra forma de extraerlos es llamando con acceso $
en el objeto a los residuos. mod1$residuals es equivalente
a residuals(mod1). Sin embargo, utilizar la función
residuals() es generalmente preferible porque es más
compatible con diferentes tipos de modelos y objetos en R.
Añadimos al gráfico el elemento de los valores estimados de y, \(\hat{y}_i\), en rojo y muestro los residuos \(\hat{\varepsilon}_i\):
ggplot(datos2, aes(x = educ, y = wage)) +
geom_point() +
geom_segment(aes(xend = educ, yend = predicciones), color = 'red', linetype = 'dashed') +
geom_point(aes(y = predicciones), color = 'red') +
geom_smooth(method = "lm", se = FALSE, color = "lightgrey") +
theme_light() +
labs(title = "Salario vs Educación con Residuos", x = "Educación", y = "Salario")## `geom_smooth()` using formula = 'y ~ x'
Podemos realizar un gráfico de dispersión para inspeccionar de forma gráfica los residuos.
plot(datos2$residuos)
El comando plot(datos2$residuos) simplemente grafica los
residuos en función de su índice, lo cual puede no ser muy informativo
(aunque ciertamente se aprecia algún patrón aglomerándose mucho por
debajo de cero, y dispersándose hacia valores más altos). Sería más útil
graficar los residuos contra los valores predichos o contra la variable
independiente para detectar patrones que indiquen violaciones de los
supuestos del modelo.
ggplot(datos2, aes(x = predicciones, y = residuos)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
theme_minimal() +
labs(title = "Residuos vs Valores Predichos", x = "Valores Predichos", y = "Residuos")
Hasta ahora vemos problemas en la dispersión de los residuos, pues no parecen estar dispersos como ruido alrededor de cero. Podríamos mejorar esto algo al ver la asimetría (skewness) de la distribución
library(e1071) # para la función skewness##
## Attaching package: 'e1071'## The following object is masked from 'package:ggplot2':
##
## elementpar(mfrow = c(1, 2)) # divide el área de gráficos en 2 columnas
plot(density(datos2$wage), main = "Gráfico de densidad: salario", ylab = "Frecuencia", sub = paste("Asimetría:", round(e1071::skewness(datos2$wage), 2))) # density para 'salario'
polygon(density(datos2$wage), col = "red")
plot(density(datos2$residuos), main = "Gráfico de densidad: residuos", ylab = "Frecuencia", sub = paste("Asimetría:", round(e1071::skewness(datos2$residuos), 2))) # density para los 'residuos'
polygon(density(datos2$residuos), col = "red")
Volveremos más tarde con la resolución de este problema, pero por ahora continuaremos con este modelo tal cual para mostrar otros elementos útiles.
4.2 Predicciones con estimados
Generaremos la predicción puntual y el intervalo de confianza para una educación promedio en años.
mean(datos2$educ)## [1] 12.56274Primero, calculamos el salario esperado para una persona con la educación promedio que es 12.56 años: \(\hat{y} = \hat{\beta}_0 + \hat{\beta}_1 \bar{x}\).
sal_pred <- ahat + bhat * 12.56274
sal_pred## [1] 5.896104La función predict() se puede igual usar para obtener la predicción puntual:
nuevo <- data.frame(educ = mean(datos2$educ))
sal_pred <- predict(mod1, newdata = nuevo)
sal_pred## 1
## 5.896103Ahora, calculamos la predicción y los intervalos para una persona con
15 años de educación para ilustrar cómo varían estos valores con
diferentes niveles de educación. Entonces (\(𝑥=15\)). El argumento
interval = prediction devuelve el valor para la predicción
puntual junto a su intervalo de predicción, que estima el rango en el
cual caerá una nueva observación individual con un cierto nivel de
confianza.
nuevo <- data.frame(educ = 15)
future_y <- predict(object = mod1, newdata = nuevo, interval = "prediction", level = 0.95)
future_y## fit lwr upr
## 1 7.215537 0.5674896 13.86358Si desean obtener solo la predicción puntual, pueden omitir el
argumento interval.
Luego, generamos el intervalo de confianza para \(𝐸(𝑦\mid𝑥)\). en este caso, debemos cambiar
el argumento a interval = "confidence". El intervalo de
confianza es para la media esperada de la variable dependiente dado un
valor específico de la independiente, aquí 15 años de educación.
future_esp_y <- predict(object = mod1, newdata = nuevo, interval = "confidence", level = 0.95)
future_esp_y <- as.data.frame(future_esp_y)
IC_inf_esp_y <- future_esp_y$lwr
IC_sup_esp_y <- future_esp_y$uprAgregando los pedazos
# Calcular intervalos de predicción para todas las observaciones
future_y_all <- predict(object = mod1, newdata = datos2, interval = "prediction", level = 0.95)
future_y_all <- as.data.frame(future_y_all)
# Calcular intervalos de confianza para todas las observaciones
future_esp_y_all <- predict(object = mod1, newdata = datos2, interval = "confidence", level = 0.95)
future_esp_y_all <- as.data.frame(future_esp_y_all)
# Combinar todos los datos en un solo data frame
nuevos_datos <- cbind(datos2, future_y_all,
IC_inf_esp_y = future_esp_y_all$lwr,
IC_sup_esp_y = future_esp_y_all$upr)Finalmente, generamos los gráficos correspondientes con los
intervalos de confianza para la predicción puntual (IC_y) y
para el valor esperado (IC_esp_y) con el siguiente
código:
IC_y <- ggplot(nuevos_datos, aes(x = educ, y = wage)) +
geom_point() +
geom_line(aes(y = lwr), color = "red", linetype = "dashed") +
geom_line(aes(y = upr), color = "red", linetype = "dashed") +
geom_smooth(method = lm, formula = y ~ x, se = TRUE, level = 0.95, col = 'blue', fill = 'pink2') +
theme_light() +
ggtitle("Predicción de y al 95%")
IC_esp_y <- ggplot(nuevos_datos, aes(x = educ, y = wage)) +
geom_point() +
geom_line(aes(y = IC_inf_esp_y), color = "blue", linetype = "dashed") +
geom_line(aes(y = IC_sup_esp_y), color = "blue", linetype = "dashed") +
geom_smooth(method = lm, formula = y ~ x, se = FALSE, col = 'blue') +
theme_light() +
ggtitle("Intervalo de Confianza de E(y|x) al 95%")Imprimimos los gráficos uno al lado del otro, para poder compararlos mejor. ¿Cuál de los dos tiene mayor amplitud?
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinegrid.arrange(IC_esp_y, IC_y, ncol = 2, nrow = 1)
Al comparar ambos gráficos, observamos que el intervalo de predicción es más amplio que el intervalo de confianza. Esto es esperado, ya que el intervalo de predicción considera la variabilidad adicional de las observaciones individuales, mientras que el intervalo de confianza se centra únicamente en la precisión de la estimación del valor medio esperado.
Veamos, a través de la comparación de dos gráficos el impacto que tiene el nivel de confianza en la amplitud de los intervalos. Para ello, tendremos que descargar e instalar algunas librerías. Se presenta primero el código y luego los gráficos obtenidos.
# Instalar y cargar las librerías necesarias si no las tiene
#install.packages("jtools")
#install.packages("gridExtra")
library(jtools)
library(gridExtra)
# Crear los gráficos con diferentes niveles de confianza
a <- plot_summs(mod1, scale = TRUE, plot.distributions = FALSE, ci_level = 0.99) +
ggtitle("Intervalos de confianza al 99%")
b <- plot_summs(mod1, scale = TRUE, plot.distributions = FALSE, ci_level = 0.95) +
ggtitle("Intervalos de confianza al 95%")
c <- plot_summs(mod1, scale = TRUE, plot.distributions = FALSE, ci_level = 0.90) +
ggtitle("Intervalos de confianza al 90%")
d <- plot_summs(mod1, scale = TRUE, plot.distributions = FALSE, ci_level = 0.70) +
ggtitle("Intervalos de confianza al 70%")
# Mostrar los gráficos uno debajo del otro
grid.arrange(a, b, c,d, ncol = 2, nrow = 2)
Al comparar los gráficos, observamos que un nivel de confianza más alto (99%) produce intervalos más amplios, ya que estamos buscando abarcar una mayor proporción de posibles valores verdaderos de los coeficientes. Por el contrario, ir bajando hacia un nivel de confianza menor (70%) resulta en intervalos más estrechos.
5 Regresión múltiple
En esta sección, queremos evaluar cómo varias variables independientes se relacionan con una variable dependiente. Específicamente, analizaremos cómo la educación y la antigüedad influyen en el salario. Luego añadiremos otras variables adicionales, como señaláramos más temprano:
wage: salario promedio por hora.
educ: años de educación.
exper: años de experiencia potencial.
tenure: años con el empleador actual (antigüedad).
nonwhite: es igual a 1 si la persona no es blanca, 0 si no.
female: es igual a 1 si la persona es mujer, 0 si no
married: es igual a 1 si la persona es casada, 0 si no.
Digamos que queremos evaluar la relación de varias variables con la dependiente. Podemos calcular las correlaciones de los pares de variables, indicando que queremos trabajar con 3 decimales:
round(cor(base1, method = "pearson"), 3)## wage educ exper antigüedad nonwhite female married lwage
## wage 1.000 0.406 0.113 0.347 -0.039 -0.340 0.229 0.937
## educ 0.406 1.000 -0.300 -0.056 -0.085 -0.085 0.069 0.431
## exper 0.113 -0.300 1.000 0.499 0.014 -0.042 0.317 0.111
## antigüedad 0.347 -0.056 0.499 1.000 0.012 -0.198 0.240 0.326
## nonwhite -0.039 -0.085 0.014 0.012 1.000 -0.011 -0.062 -0.039
## female -0.340 -0.085 -0.042 -0.198 -0.011 1.000 -0.166 -0.374
## married 0.229 0.069 0.317 0.240 -0.062 -0.166 1.000 0.271
## lwage 0.937 0.431 0.111 0.326 -0.039 -0.374 0.271 1.000
## expersq 0.030 -0.331 0.961 0.459 0.009 -0.028 0.217 0.023
## antigüedadcuad 0.267 -0.069 0.423 0.922 -0.007 -0.176 0.167 0.236
## expersq antigüedadcuad
## wage 0.030 0.267
## educ -0.331 -0.069
## exper 0.961 0.423
## antigüedad 0.459 0.922
## nonwhite 0.009 -0.007
## female -0.028 -0.176
## married 0.217 0.167
## lwage 0.023 0.236
## expersq 1.000 0.414
## antigüedadcuad 0.414 1.000Antes de estimar el modelo de regresión múltiple, es útil explorar las relaciones entre las variables mediante la matriz de correlaciones. Esto nos ayuda a detectar posibles problemas de multicolinealidad y a entender las relaciones bivariadas. En este caso, haremos una restricción a tres variables, ayudando a mejorar la legibilidad.
# Seleccionamos las variables de interés
variables_interes <- base1[, c("wage", "educ", "antigüedad")]
# Calculamos la matriz de correlaciones y redondeamos a 3 decimales
matriz_correlaciones <- round(cor(variables_interes, method = "pearson"), 3)
matriz_correlaciones## wage educ antigüedad
## wage 1.000 0.406 0.347
## educ 0.406 1.000 -0.056
## antigüedad 0.347 -0.056 1.000Interpretación:
- Salario y Educación: Un coeficiente positivo indica que a mayor nivel educativo, el salario tiende a ser mayor. La fuerza de esta relación parece ser moderada, 0.406.
- Salario y Antigüedad: Un coeficiente positivo sugiere que con más años de antigüedad, el salario también aumenta. Parece ser de una fuerza similar a educación.
- Educación y Antigüedad: El coeficiente es negativo, y pequeño. Entre estas variables la relación no es clara en fuerza pero parecen ir en direcciones opuestas. Es menos probable que haya multicolinealidad entre estas variables independientes.
Para visualizar la relación conjunta entre las tres variables, utilizamos gráficos tridimensionales.
library(ggplot2)
library(plotly)
attach(base1)
plot_ly(x = educ, y = antigüedad, z = wage, type = "scatter3d", color = wage) |>
layout(scene = list(xaxis = list(title = 'educación (en años)'),
yaxis = list(title = 'antigüedad (en años)'),
zaxis = list(title = 'Salario (en USD/h)')))## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-modeEste gráfico interactivo permite rotar y explorar la relación entre las variables desde diferentes ángulos, facilitando la identificación de patrones. El próximo es menos interactivo, pero logra un efecto similar:
library(scatterplot3d)
graf <- scatterplot3d(x = educ, y = antigüedad, z = wage, pch = 16,
cex.lab = 1, highlight.3d = TRUE, type = "h",
xlab = 'Años de educación',
ylab = 'Antigüedad (años)',
zlab = 'Salario (USD/h)')
En este caso, las líneas verticales ayudan a visualizar la posición de cada punto en el espacio tridimensional.
Ahora, estimamos un modelo de regresión lineal múltiple para cuantificar el efecto de la educación y la antigüedad en el salario.
# Estimamos el modelo de regresión múltiple
mod2 <- lm(wage ~ educ + antigüedad, data = base1)
# Resumen del modelo
summary(mod2)##
## Call:
## lm(formula = wage ~ educ + antigüedad, data = base1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.1438 -1.7288 -0.6372 1.2575 14.7482
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.22162 0.64015 -3.47 0.000563 ***
## educ 0.56914 0.04881 11.66 < 2e-16 ***
## antigüedad 0.18958 0.01871 10.13 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.092 on 523 degrees of freedom
## Multiple R-squared: 0.3019, Adjusted R-squared: 0.2992
## F-statistic: 113.1 on 2 and 523 DF, p-value: < 2.2e-16La salida es similar al modelo anterior, con filas adicionales para cada coeficiente estimado adicional. El intercepto o constante representa el salario promedio cuando la educación y la antigüedad son cero (interpretación limitada en este contexto). Educ indica el cambio esperado en el salario por cada año adicional de educación, manteniendo constante la antigüedad en su valor esperado. Antigüedad muestra el cambio esperado en el salario por cada año adicional de antigüedad, manteniendo constante la educación. En este caso vemos que un año adicional de educación se traduce a unos 57 chavos adicionales en ingreso, mientras que antigüedad añade cerca de 19 chavos al ingreso.
Todos los coeficientes reportan valores t relativamente grandes en relación a sus errores estándar, y rechazamos en cada caso la hipótesis nula de no ser significativamente distinto a un efecto nulo. Notamos que tanto \(R^2\) y \(R^2\) ajustado han aumentado ambos, aunque la distancia entre ambos ahora es algo más notable: duplicamos la varianza explicada.
¿Cómo se ve esto? Añadimos el plano de regresión al gráfico 3D para visualizar cómo el modelo ajusta los datos.
graf <- scatterplot3d(x = educ, y = antigüedad, z = wage, pch = 16,
cex.lab = 1, highlight.3d = TRUE, type = "h",
xlab = 'Años de educación',
ylab = 'Experiencia (años)',
zlab = 'Salario (USD/h)')
graf$plane3d(mod2, lty.box = "solid", col = 'mediumblue')
El plano representa las predicciones del modelo para diferentes combinaciones de educación y antigüedad. Podemos observar qué tan bien el plano ajusta a los datos reales, intentando cortar por un espacio vectorial estimado linealmente.
El ANOVA nos permite evaluar la significancia global del modelo y la contribución de cada variable independiente.
anova(mod2)## Analysis of Variance Table
##
## Response: wage
## Df Sum Sq Mean Sq F value Pr(>F)
## educ 1 1179.7 1179.73 123.43 < 2.2e-16 ***
## antigüedad 1 981.7 981.72 102.71 < 2.2e-16 ***
## Residuals 523 4999.0 9.56
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Sumas de Cuadrados (Sum Sq): Indican la variabilidad explicada por cada variable independiente y la variabilidad residual.
Para garantizar la validez de nuestro modelo de regresión múltiple, es fundamental verificar que se cumplen los supuestos básicos del modelo lineal. A continuación, revisaremos cada uno de estos supuestos en el orden mencionado anteriormente, aplicando pruebas diagnósticas y proporcionando interpretaciones detalladas.
5.1 Análisis de suposiciones
Antes evaluamos elementos de los residuos, pero lo hicimos por encima. R tiene en forma base unos gráficos que ayudan a informarnos sobre si los modelos tienen problemas en cómo están siendo aplicados.
Primero, podemos evaluar el modelo gráficamente con la función
plot():
plot(mod2)
- Residuos vs Ajustados (Residuals vs Fitted):
Qué muestra: Este gráfico representa los residuos estandarizados en función de los valores ajustados por el modelo.
Interpretación: Sirve para detectar patrones no lineales y evaluar la homogeneidad de la varianza (homocedasticidad). Si los puntos se distribuyen aleatoriamente alrededor de la línea horizontal (residuo = 0) sin formar patrones, indica que el modelo es adecuado. Patrones sistemáticos o formas específicas (como una curva) sugieren que el modelo no captura adecuadamente la relación entre las variables, o que existe heterocedasticidad.
- Gráfico Q-Q Normal (Normal Q-Q Plot):
- Qué muestra: Compara la distribución de los residuos estandarizados con una distribución normal teórica.
- Interpretación: Evalúa la normalidad de los residuos, un supuesto clave en modelos lineales. Si los puntos siguen aproximadamente una línea recta, los residuos se distribuyen normalmente. Desviaciones significativas de la línea recta indican que los residuos no son normales, lo que puede afectar la validez de los intervalos de confianza y pruebas de hipótesis.
- Escala-Ubicación (Scale-Location Plot):
- Qué muestra: Grafica la raíz cuadrada de los residuos estandarizados
(sqrt(\|Residuos estandarizados\|))frente a los valores ajustados. - Interpretación: Ayuda a verificar la homocedasticidad. Una dispersión uniforme de puntos sugiere varianza constante de los residuos. Si los puntos muestran un patrón (por ejemplo, se ensanchan o estrechan a lo largo del eje de los ajustados), indica heterocedasticidad, lo que puede afectar la eficiencia de los estimadores.
- Residuos Estandarizados vs Apalancamiento (Residuals vs Leverage):
- Qué muestra: Muestra los residuos estandarizados frente al apalancamiento de cada observación, con curvas de distancia de Cook superpuestas.
- Interpretación: Identifica observaciones influyentes que tienen un gran impacto en el ajuste del modelo. Puntos con alto apalancamiento y residuos grandes pueden distorsionar los resultados. Las líneas de distancia de Cook ayudan a detectar estos puntos. Observaciones más allá de estas líneas merecen una revisión adicional.
Notamos que el modelo tiene sus problemas, violentando varios de los principios señalados antes.
5.1.1 Linealidad
La relación entre las variables independientes y la variable dependiente es lineal.
# Gráficos de dispersión con línea de regresión para cada variable independiente
library(ggplot2)
ggplot(base1, aes(x = educ, y = wage)) +
geom_point() +
geom_smooth(method = "lm", col = "red") +
labs(title = "Salario vs Educación", x = "Años de Educación", y = "Salario") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
ggplot(base1, aes(x = antigüedad, y = wage)) +
geom_point() +
geom_smooth(method = "lm", col = "red") +
labs(title = "Salario vs Antigüedad", x = "Años de Antigüedad", y = "Salario") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
5.1.2 Independencia de errores
Realizando la prueba de Durbin-Watson para verificar la autocorrelación estocástica. La hipótesis nula es que los residuos no tienen autocorrelación (es posible una alternativa, editando opciones), es decir que los residuos son normales.
# Prueba de Durbin-Watson
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericdwtest(mod2)##
## Durbin-Watson test
##
## data: mod2
## DW = 1.7907, p-value = 0.008009
## alternative hypothesis: true autocorrelation is greater than 0En este caso la prueba rechaza la nula. Tenemos problemas en nuestros residuos tal cual modelados al presente, pues se viola el supuesto de independencia.
5.1.3 Homoscedasticidad
# Gráfico de Residuos vs Valores Ajustados
ggplot(data = base1, aes(x = mod2$fitted.values, y = mod2$residuals)) +
geom_point() +
geom_smooth(method = "loess", col = "red", se = FALSE) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "Residuos vs Valores Ajustados", x = "Valores Ajustados", y = "Residuos") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
# Prueba de Breusch-Pagan
library(lmtest)
bptest(mod2)##
## studentized Breusch-Pagan test
##
## data: mod2
## BP = 40.798, df = 2, p-value = 1.383e-09Vemos una representación visual (con ggplot()) de los
residuos contrastados a los valores ajustados. La prueba Breusch Pagan
tiene la hipótesis nula siguiente: varianza constante. Por consiguiente,
rechazarla es encontrar heteroscedasticidad.
5.1.4 Normalidad de errores
El supuesto a verificar ahora (aunque no está yendo bien en general) es saber si los errores se distribuyen normalmente. Hay varias opciones gráficas y con pruebas estadísticas: el gráfico Q-Q Plot, el histograma de los residuos, el diagrama de densidad de residuos, y la prueba Shapiro-Wilk.
# Gráfico Q-Q de los residuos
qqnorm(mod2$residuals)
qqline(mod2$residuals, col = "red")
# Histograma de los residuos
hist(mod2$residuals, breaks = 20, main = "Histograma de Residuos", xlab = "Residuos", col = "lightblue")
# Gráfico de densidad de los residuos con asimetría
library(e1071) # Para la función skewness()
plot(density(mod2$residuals),
main = "Gráfico de Densidad de Residuos", xlab = "Residuos",
sub = paste("Asimetría:", round(skewness(mod2$residuals), 2)))
polygon(density(mod2$residuals), col = "red")
Por ejemplo, estos residuos distan de ser normales. Las cuantilas de los residuos distan de la normalidad esperada, por MUCHO. Verificamos
shapiro.test(mod2$residuals)##
## Shapiro-Wilk normality test
##
## data: mod2$residuals
## W = 0.88987, p-value < 2.2e-16La hipótesis nula de la prueba de Shapiro-Wilk es que estamos ante normalidad en los errores. Esto confirma lo sugerido por el análisis gráfico y estadístico precedente.
En los siguientes gráficos se muestran los residuos contra cada uno de los regresores, los cuales se realizan con el siguiente código:
plot1 <- ggplot(data = base1, aes(educ, mod1$residuals)) +
geom_point() +
geom_smooth(color = "firebrick") +
geom_hline(yintercept = 0) +
theme_bw()
plot2 <- ggplot(data = base1, aes(antigüedad, mod1$residuals)) +
geom_point() +
geom_smooth(color = "firebrick") +
geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot1, plot2)## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
El análisis gráfico indica problemas de heteroscedasticidad y de correlación entre los residuos y el nivel de los regresores.
5.1.5 Multicolinealidad
# Matriz de correlaciones entre variables independientes
vars_indep <- base1[, c("educ", "antigüedad")]
cor(vars_indep)## educ antigüedad
## educ 1.00000000 -0.05617257
## antigüedad -0.05617257 1.00000000# Cálculo del VIF
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## somevif(mod2)## educ antigüedad
## 1.003165 1.003165Los valores no son particularmente altos en correlación (alto es cerca de 1 ó -1, bajo es cerca de 0). En nuestro caso, la correlación entre educación y antigüedad es baja (-0.056), lo que sugiere baja multicolinealidad. La prueba del factor de inflación de varianza (VIF) se puede usar para verificar numéricamente si hay multicolinealidad.
5.1.6 Identificar observaciones influyentes
Aunque no es un supuesto básico, es importante identificar puntos que puedan influir excesivamente en el modelo. Esto se puede verificar: - Analizando el gráfico de Residuos Estandarizados vs Apalancamiento. - Calculando la Distancia de Cook.
# Gráfico de Residuos Estandarizados vs Apalancamiento
plot(mod2, which = 5)
# Gráfico de Distancia de Cook
plot(mod2, which = 4)
# Identificar observaciones con alta Distancia de Cook
cooksd <- cooks.distance(mod2)
influential <- as.numeric(names(cooksd)[(cooksd > (4 / nrow(base1)))])
influential## [1] 10 15 16 58 59 66 74 98 112 128 170 171 186 203 229 237 245 260 265
## [20] 278 282 297 307 326 343 379 398 406 440 458 463 464 465 476 498 503 505 522
## [39] 524# Mostrar las observaciones influyentes
base1[influential, ]## wage educ exper antigüedad nonwhite female married lwage expersq
## 10 18.18 17 22 21 0 0 1 2.9003222 484
## 15 22.20 12 31 15 0 0 1 3.1000924 961
## 16 17.33 16 14 0 0 0 1 2.8524392 196
## 58 10.00 8 13 0 1 0 0 2.3025851 169
## 59 21.63 18 8 8 0 1 0 3.0740812 64
## 66 19.98 14 26 23 0 0 1 2.9947317 676
## 74 2.91 12 20 34 0 1 1 1.0681531 400
## 98 13.16 12 34 22 0 0 1 2.5771818 1156
## 112 24.98 18 29 25 0 0 1 3.2180755 841
## 128 1.50 8 31 30 0 0 0 0.4054651 961
## 170 15.38 13 25 21 0 0 1 2.7330680 625
## 171 14.58 18 13 7 0 1 0 2.6796508 169
## 186 21.86 12 24 16 0 0 1 3.0846586 576
## 203 10.00 8 9 0 0 0 1 2.3025851 81
## 229 22.86 16 16 7 0 0 1 3.1293886 256
## 237 3.50 18 3 1 0 1 0 1.2527629 9
## 245 18.16 16 29 7 0 0 1 2.8992214 841
## 260 18.00 18 13 0 0 1 0 2.8903718 169
## 265 8.75 12 47 44 0 0 1 2.1690538 2209
## 278 18.89 17 26 20 0 0 1 2.9386327 676
## 282 2.95 14 41 23 1 0 1 1.0818052 1681
## 297 3.60 10 34 26 0 1 1 1.2809339 1156
## 307 6.25 11 35 31 0 0 1 1.8325815 1225
## 326 18.16 16 35 28 0 0 1 2.8992214 1225
## 343 15.00 11 35 31 0 0 1 2.7080503 1225
## 379 4.17 0 22 10 0 1 0 1.4279160 484
## 398 4.29 12 47 25 0 1 1 1.4562867 2209
## 406 8.43 8 27 3 0 0 1 2.1317968 729
## 440 20.00 12 22 4 0 0 1 2.9957323 484
## 458 6.50 14 41 33 0 0 1 1.8718022 1681
## 463 3.05 8 50 24 0 1 0 1.1151416 2500
## 464 3.50 12 26 20 0 1 1 1.2527629 676
## 465 2.92 3 51 30 1 0 0 1.0715836 2601
## 476 17.50 16 23 22 0 0 1 2.8622010 529
## 498 5.15 16 39 21 0 1 0 1.6389967 1521
## 503 2.89 0 42 0 0 1 1 1.0612565 1764
## 505 17.71 16 10 3 0 0 1 2.8741293 100
## 522 15.00 16 14 2 0 1 1 2.7080503 196
## 524 4.67 15 13 18 0 0 1 1.5411590 169
## antigüedadcuad
## 10 441
## 15 225
## 16 0
## 58 0
## 59 64
## 66 529
## 74 1156
## 98 484
## 112 625
## 128 900
## 170 441
## 171 49
## 186 256
## 203 0
## 229 49
## 237 1
## 245 49
## 260 0
## 265 1936
## 278 400
## 282 529
## 297 676
## 307 961
## 326 784
## 343 961
## 379 100
## 398 625
## 406 9
## 440 16
## 458 1089
## 463 576
## 464 400
## 465 900
## 476 484
## 498 441
## 503 0
## 505 9
## 522 4
## 524 324Puntos fuera de las líneas de referencia pueden ser influyentes. Los gráficos señalan los valores que parecen distantes al resto de lo observado. En estos casos, se recomienda revisar las observaciones para verificar si hay errores en los datos o si representan casos especiales.
5.1.7 Análisis de residuos contra variables independientes
Es útil graficar los residuos contra cada una de las variables independientes para detectar patrones específicos.
# Residuos vs Educación
plot1 <- ggplot(data = base1, aes(x = educ, y = mod2$residuals)) +
geom_point() +
geom_smooth(method = "loess", color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "Residuos vs Educación", x = "Educación", y = "Residuos") +
theme_bw()
# Residuos vs Antigüedad
plot2 <- ggplot(data = base1, aes(x = antigüedad, y = mod2$residuals)) +
geom_point() +
geom_smooth(method = "loess", color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "Residuos vs Antigüedad", x = "Antigüedad", y = "Residuos") +
theme_bw()
# Mostrar los gráficos lado a lado
library(gridExtra)
grid.arrange(plot1, plot2, ncol = 2)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
- Residuos vs Educación: Si observamos un patrón sistemático, podría indicar que la relación no es completamente lineal o que hay variables omitidas.
- Residuos vs Antigüedad: Patrones similares pueden sugerir problemas de heterocedasticidad o relaciones no lineales.
5.1.8 ¿Qué hacemos?
Como detectamos problemas de heterocedasticidad y normalidad de errores, podemos probar transformando la variable dependiente. Por ahora añadiré unas variables adicionales y también haremos unas variaciones con transformaciones (una de estas ya estaba en el conjunto de datos). Es posible que al incorporar variables relevantes que puedan explicar mejor la variabilidad en el salario, mejore el ajuste del modelo. Por otro lado, es posible que transformar la variable dependiente o algunas independientes para corregir violaciones a los supuestos del modelo lineal arreglen estos problemas.
A continuación, implementaremos estas estrategias paso a paso.
Añadiremos variables que podrían influir significativamente en el salario, como:
- married: Estado civil (1 si está casado, 0 si no).
- female: Género (1 si es mujer, 0 si es hombre).
- nonwhite: Etnicidad (1 si no es blanco, 0 si es blanco).
# Actualizamos el modelo añadiendo la variable 'married'
mod3 <- update(mod2, wage ~ educ + antigüedad + married)
# Resumen del modelo
summary(mod3)##
## Call:
## lm(formula = wage ~ educ + antigüedad + married, data = base1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.0742 -1.7278 -0.5579 1.1913 14.5374
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.52575 0.64205 -3.934 9.49e-05 ***
## educ 0.55614 0.04857 11.450 < 2e-16 ***
## antigüedad 0.17482 0.01913 9.139 < 2e-16 ***
## married 0.89235 0.28310 3.152 0.00171 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.066 on 522 degrees of freedom
## Multiple R-squared: 0.3149, Adjusted R-squared: 0.311
## F-statistic: 79.98 on 3 and 522 DF, p-value: < 2.2e-16- Intercepto (-2.52575): Representa el salario promedio cuando educ, antigüedad y married son cero. Aunque no es interpretable en este contexto, es necesario para el modelo.
- Educación (0.55614): Por cada año adicional de educación, el salario promedio aumenta en aproximadamente 56 chavos, manteniendo constantes la antigüedad y el estado civil.
- Antigüedad (0.17482): Por cada año adicional de antigüedad con el empleador actual, el salario promedio aumenta en aproximadamente $0.17, manteniendo constantes las otras variables.
- Casado (0.89235): Estar casado se asocia con un aumento promedio de 89 chavos en el salario, manteniendo constantes la educación y la antigüedad.
- Significancia estadística: Todos los coeficientes son estadísticamente significativos (p < 0.05), lo que indica que tienen un efecto significativo en el salario.
- Ajuste del modelo:
- R-cuadrado (0.3149): El modelo explica aproximadamente el 31% de la variabilidad en el salario, lo que es una mejora respecto a los modelos anteriores.
- Error estándar residual (3.066): Es menor que en los modelos previos, lo que indica un mejor ajuste.
Para visualizar los coeficientes del modelo, podemos utilizar la
librería modelsummary y la función
modelplot():
library(modelsummary)
modelplot(mod3, coef_omit = "Intercept", color = "blue", size = 1) +
labs(title = "Coeficientes del Modelo 3")
Nota: Al incluir variables categóricas como married, es importante considerar que sus coeficientes representan diferencias respecto a la categoría de referencia.
5.1.8.1 Estandarización de Variables
Para comparar los efectos de las variables en la misma escala, podemos estandarizar las variables numéricas:
base1_estandarizado <- base1 |>
mutate(across(where(is.numeric), scale))
# Ajustamos el modelo con los datos estandarizados
mod3_est <- lm(wage ~ educ + antigüedad + married, data = base1_estandarizado)
# Visualicemos los coeficientes estandarizados
modelplot(mod3_est, coef_omit = "Intercept", color = "blue", size = 1) +
labs(title = "Coeficientes Estandarizados del Modelo 3")
Los coeficientes estandarizados permiten comparar directamente el efecto relativo de cada variable en el salario. Un coeficiente más grande en valor absoluto indica un efecto mayor.
Al ejecutar plot(mod3), obtenemos los gráficos
diagnósticos para evaluar los supuestos del modelo.
plot(mod3)
#library(lmtest)
bptest(mod3)##
## studentized Breusch-Pagan test
##
## data: mod3
## BP = 43.296, df = 3, p-value = 2.129e-09shapiro.test(mod3$residuals)##
## Shapiro-Wilk normality test
##
## data: mod3$residuals
## W = 0.89096, p-value < 2.2e-16Los residuos todavía siguen distando de ser una distribución normal. Ampliaremos el modelo con otras variables que señaláramos antes como de interés.
# Creamos un nuevo modelo incluyendo 'female' y 'nonwhite'
mod4 <- lm(wage ~ educ + antigüedad + married + female + nonwhite, data = base1)
# Resumen del modelo
summary(mod4)##
## Call:
## lm(formula = wage ~ educ + antigüedad + married + female + nonwhite,
## data = base1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.4948 -1.7690 -0.5398 1.0261 13.9814
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.12270 0.66394 -1.691 0.0914 .
## educ 0.52878 0.04718 11.208 < 2e-16 ***
## antigüedad 0.15423 0.01873 8.234 1.47e-15 ***
## married 0.68258 0.27546 2.478 0.0135 *
## female -1.71149 0.26644 -6.424 3.01e-10 ***
## nonwhite -0.06516 0.42715 -0.153 0.8788
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.956 on 520 degrees of freedom
## Multiple R-squared: 0.3653, Adjusted R-squared: 0.3592
## F-statistic: 59.85 on 5 and 520 DF, p-value: < 2.2e-16- Educación y antigüedad: Siguen siendo significativas y positivas.
- Married (0.68258): Estar casado se asocia con un aumento promedio de $0.68 en el salario, manteniendo constantes las demás variables.
- Female (-1.71149): Ser mujer se asocia con una disminución promedio de $1.71 en el salario, manteniendo constantes las demás variables. Este efecto es estadísticamente significativo.
- Nonwhite (-0.06516): No es estadísticamente significativo (p = 0.8788), lo que sugiere que, en este modelo, la variable nonwhite no tiene un efecto significativo en el salario.
- Ajuste del modelo:
- R-cuadrado (0.3653), ajustado (0.3592): El modelo explica aproximadamente el 36% de la variabilidad en el salario, mejorando levemente respecto al modelo anterior.
5.1.8.2 Análisis
Es importante verificar la multicolinealidad al añadir nuevas variables.
#library(car)
vif(mod4)## educ antigüedad married female nonwhite
## 1.025155 1.099796 1.087978 1.066217 1.011546No parece haber problema de multicolinealidad, todas están cerca de 1.
Verificaremos rápido igual los gráficos del modelo.
plot(mod4)
#library(lmtest)
bptest(mod4)##
## studentized Breusch-Pagan test
##
## data: mod4
## BP = 43.589, df = 5, p-value = 2.806e-08shapiro.test(mod4$residuals)##
## Shapiro-Wilk normality test
##
## data: mod4$residuals
## W = 0.88709, p-value < 2.2e-16Seguimos rechazando la hipótesis nula de Shapiro-Wilk aún en este caso. El problema no parece ser resuelto por añadidura de otros términos. Aplicamos una transformación logarítmica a la variable dependiente, salario.
# Estimamos el modelo con el logaritmo del salario
mod5 <- lm(log(wage) ~ educ + antigüedad + married + female + nonwhite, data = base1)
# Resumen del modelo
summary(mod5)##
## Call:
## lm(formula = log(wage) ~ educ + antigüedad + married + female +
## nonwhite, data = base1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.89855 -0.27105 -0.03577 0.24487 1.24663
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.570835 0.092937 6.142 1.62e-09 ***
## educ 0.079483 0.006604 12.035 < 2e-16 ***
## antigüedad 0.019443 0.002622 7.416 4.94e-13 ***
## married 0.146703 0.038558 3.805 0.000159 ***
## female -0.280401 0.037296 -7.518 2.45e-13 ***
## nonwhite -0.002419 0.059791 -0.040 0.967738
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4138 on 520 degrees of freedom
## Multiple R-squared: 0.3996, Adjusted R-squared: 0.3939
## F-statistic: 69.23 on 5 and 520 DF, p-value: < 2.2e-16Interpretación:
Coeficientes interpretados en términos logarítmicos: - Educación (0.079483): Un año adicional de educación se asocia con un aumento promedio del 7.95% en el salario, manteniendo constantes las demás variables. - Antigüedad (0.019443): Un año adicional de antigüedad se asocia con un aumento promedio del 1.94% en el salario. - Married (0.146703): Estar casado se asocia con un aumento promedio del 14.67% en el salario. - Female (-0.280401): Ser mujer se asocia con una disminución promedio del 28.04% en el salario. - Nonwhite (-0.002419): No es significativo. - Mejora en los supuestos del modelo: - La transformación logarítmica puede ayudar a corregir problemas de heterocedasticidad y normalidad de los residuos.
plot(mod5)
Estos gráficos se ven distintos a los anteriores. Los residuos parecen dispersarse de manera más uniforme alrededor de cero. Hagamos varias pruebas para verificar:
#library(lmtest)
bptest(mod5)##
## studentized Breusch-Pagan test
##
## data: mod5
## BP = 11.928, df = 5, p-value = 0.03579shapiro.test(mod5$residuals)##
## Shapiro-Wilk normality test
##
## data: mod5$residuals
## W = 0.98592, p-value = 5.75e-05Sigue rechazándose la nula, aunque por menos margen que antes en el caso de Breusch-Pagan (y menos mejorado en Shapiro-Wilk).
Haremos unos últimos intentos de ajuste, como el uso de errores estándar robustos (ajustan las estimaciones de la varianza para corregir la heterocedasticidad sin cambiar los coeficientes estimados).
# Instalar y cargar los paquetes necesarios
#install.packages("lmtest")
#install.packages("sandwich")
library(lmtest)
library(sandwich)
# Recalcular los errores estándar usando la matriz de varianza-covarianza robusta
coef(summary(mod5))## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.570834892 0.092936693 6.14219071 1.618659e-09
## educ 0.079482847 0.006604044 12.03548155 1.330233e-29
## antigüedad 0.019442804 0.002621755 7.41594863 4.940979e-13
## married 0.146703436 0.038558103 3.80473683 1.588325e-04
## female -0.280400745 0.037295927 -7.51826714 2.448824e-13
## nonwhite -0.002419444 0.059791466 -0.04046471 9.677382e-01coeftest(mod5, vcov = vcovHC(mod5, type = "HC1"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5708349 0.1019869 5.5971 3.531e-08 ***
## educ 0.0794828 0.0072800 10.9180 < 2.2e-16 ***
## antigüedad 0.0194428 0.0031956 6.0842 2.272e-09 ***
## married 0.1467034 0.0400634 3.6618 0.0002761 ***
## female -0.2804007 0.0378265 -7.4128 5.048e-13 ***
## nonwhite -0.0024194 0.0614213 -0.0394 0.9685938
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Los coeficientes estimados permanecen iguales en este caso, pero los errores estándar y los valores p pueden cambiar. Para más información sobre este método recomiendo este artículo.
Si quisiéramos hacer otras transformaciones adicionales, no satisfechos con estos errores robustos, podríamos buscar una transformación Box-Cox, un tipo de transformación de potencia que redistribuye la variable con un logaritmo y luego eleva a un exponente óptimo para normalizar la variable. Para más información recomiendo este enlace.
# Instalar y cargar el paquete MASS
#install.packages("MASS")
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:plotly':
##
## select## The following object is masked from 'package:dplyr':
##
## select## The following object is masked from 'package:wooldridge':
##
## cement# Encontrar el lambda óptimo para la transformación Box-Cox
boxcox_mod <- boxcox(mod4, plotit = TRUE)
lambda_optimo <- boxcox_mod$x[which.max(boxcox_mod$y)]
lambda_optimo## [1] -0.1414141# Aplicar la transformación Box-Cox al salario
base1$wage_boxcox <- (base1$wage^lambda_optimo - 1) / lambda_optimo
# Reestimar el modelo con la variable transformada
mod_boxcox <- lm(wage_boxcox ~ educ + antigüedad + married + female + nonwhite, data = base1)
summary(mod_boxcox)##
## Call:
## lm(formula = wage_boxcox ~ educ + antigüedad + married + female +
## nonwhite, data = base1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.81512 -0.21056 -0.02169 0.20059 0.90536
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.612834 0.073483 8.340 6.71e-16 ***
## educ 0.062016 0.005222 11.877 < 2e-16 ***
## antigüedad 0.014772 0.002073 7.126 3.47e-12 ***
## married 0.120575 0.030487 3.955 8.72e-05 ***
## female -0.220928 0.029489 -7.492 2.94e-13 ***
## nonwhite -0.001225 0.047276 -0.026 0.979
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3272 on 520 degrees of freedom
## Multiple R-squared: 0.3941, Adjusted R-squared: 0.3883
## F-statistic: 67.64 on 5 and 520 DF, p-value: < 2.2e-16# Verificar los supuestos nuevamente
shapiro.test(mod_boxcox$residuals)##
## Shapiro-Wilk normality test
##
## data: mod_boxcox$residuals
## W = 0.98056, p-value = 1.786e-06bptest(mod_boxcox)##
## studentized Breusch-Pagan test
##
## data: mod_boxcox
## BP = 8.7255, df = 5, p-value = 0.1205No rechazamos la hipótesis nula de homocedasticidad por Breusch-Pagan. La heterocedasticidad se ha mitigado en este modelo. Sin embargo, rechazamos la hipótesis nula de normalidad. Los residuos aún no siguen una distribución normal.
También podríamos incluir otros términos polinomiales o interacciones.
# Modelo con términos cuadráticos
mod_poly <- lm(log(wage) ~ educ + I(educ^2) + antigüedad + I(antigüedad^2) + married + female + nonwhite, data = base1)
# Resumen del modelo
summary(mod_poly)##
## Call:
## lm(formula = log(wage) ~ educ + I(educ^2) + antigüedad + I(antigüedad^2) +
## married + female + nonwhite, data = base1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.86343 -0.25631 -0.02747 0.24719 1.26424
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.0890157 0.1836224 5.931 5.52e-09 ***
## educ -0.0200558 0.0298881 -0.671 0.502501
## I(educ^2) 0.0042157 0.0012474 3.380 0.000780 ***
## antigüedad 0.0407787 0.0064974 6.276 7.34e-10 ***
## I(antigüedad^2) -0.0008442 0.0002312 -3.652 0.000287 ***
## married 0.1264775 0.0381141 3.318 0.000969 ***
## female -0.2668000 0.0367257 -7.265 1.38e-12 ***
## nonwhite -0.0262610 0.0586796 -0.448 0.654678
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4048 on 518 degrees of freedom
## Multiple R-squared: 0.4277, Adjusted R-squared: 0.42
## F-statistic: 55.31 on 7 and 518 DF, p-value: < 2.2e-16# Verificar los supuestos
shapiro.test(mod_poly$residuals)##
## Shapiro-Wilk normality test
##
## data: mod_poly$residuals
## W = 0.98657, p-value = 9.058e-05bptest(mod_poly)##
## studentized Breusch-Pagan test
##
## data: mod_poly
## BP = 11.462, df = 7, p-value = 0.1197Tampoco rechazamos en esta transformación la hipótesis nula de homoscedasticidad. La heteroscedasticidad se ha reducido. Pero sí volvemos a notar que los errores no son normales. Soluciones como bootstrapping podrían ser útiles, así como modelos generalizados. Por ahora, podemos enfocarnos en sacar las tablas de los modelos, para uso en o para otros programas.
library(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerstargazer(mod1,mod2,mod3,mod4,mod5,mod_boxcox) #de base, para LaTeX##
## % Table created by stargazer v.5.2.3 by Marek Hlavac, Social Policy Institute. E-mail: marek.hlavac at gmail.com
## % Date and time: Fri, Oct 10, 2025 - 12:48:58
## \begin{table}[!htbp] \centering
## \caption{}
## \label{}
## \begin{tabular}{@{\extracolsep{5pt}}lcccccc}
## \\[-1.8ex]\hline
## \hline \\[-1.8ex]
## & \multicolumn{6}{c}{\textit{Dependent variable:}} \\
## \cline{2-7}
## \\[-1.8ex] & \multicolumn{4}{c}{wage} & log(wage) & wage\_boxcox \\
## \\[-1.8ex] & (1) & (2) & (3) & (4) & (5) & (6)\\
## \hline \\[-1.8ex]
## educ & 0.541$^{***}$ & 0.569$^{***}$ & 0.556$^{***}$ & 0.529$^{***}$ & 0.079$^{***}$ & 0.062$^{***}$ \\
## & (0.053) & (0.049) & (0.049) & (0.047) & (0.007) & (0.005) \\
## & & & & & & \\
## antigüedad & & 0.190$^{***}$ & 0.175$^{***}$ & 0.154$^{***}$ & 0.019$^{***}$ & 0.015$^{***}$ \\
## & & (0.019) & (0.019) & (0.019) & (0.003) & (0.002) \\
## & & & & & & \\
## married & & & 0.892$^{***}$ & 0.683$^{**}$ & 0.147$^{***}$ & 0.121$^{***}$ \\
## & & & (0.283) & (0.275) & (0.039) & (0.030) \\
## & & & & & & \\
## female & & & & $-$1.711$^{***}$ & $-$0.280$^{***}$ & $-$0.221$^{***}$ \\
## & & & & (0.266) & (0.037) & (0.029) \\
## & & & & & & \\
## nonwhite & & & & $-$0.065 & $-$0.002 & $-$0.001 \\
## & & & & (0.427) & (0.060) & (0.047) \\
## & & & & & & \\
## Constant & $-$0.905 & $-$2.222$^{***}$ & $-$2.526$^{***}$ & $-$1.123$^{*}$ & 0.571$^{***}$ & 0.613$^{***}$ \\
## & (0.685) & (0.640) & (0.642) & (0.664) & (0.093) & (0.073) \\
## & & & & & & \\
## \hline \\[-1.8ex]
## Observations & 526 & 526 & 526 & 526 & 526 & 526 \\
## R$^{2}$ & 0.165 & 0.302 & 0.315 & 0.365 & 0.400 & 0.394 \\
## Adjusted R$^{2}$ & 0.163 & 0.299 & 0.311 & 0.359 & 0.394 & 0.388 \\
## Residual Std. Error & 3.378 (df = 524) & 3.092 (df = 523) & 3.066 (df = 522) & 2.956 (df = 520) & 0.414 (df = 520) & 0.327 (df = 520) \\
## F Statistic & 103.363$^{***}$ (df = 1; 524) & 113.067$^{***}$ (df = 2; 523) & 79.978$^{***}$ (df = 3; 522) & 59.848$^{***}$ (df = 5; 520) & 69.227$^{***}$ (df = 5; 520) & 67.641$^{***}$ (df = 5; 520) \\
## \hline
## \hline \\[-1.8ex]
## \textit{Note:} & \multicolumn{6}{r}{$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01} \\
## \end{tabular}
## \end{table}stargazer(mod1,mod2,mod3, type="text", title = "Resultados de regresión")##
## Resultados de regresión
## =============================================================================================
## Dependent variable:
## -------------------------------------------------------------------------
## wage
## (1) (2) (3)
## ---------------------------------------------------------------------------------------------
## educ 0.541*** 0.569*** 0.556***
## (0.053) (0.049) (0.049)
##
## antigüedad 0.190*** 0.175***
## (0.019) (0.019)
##
## married 0.892***
## (0.283)
##
## Constant -0.905 -2.222*** -2.526***
## (0.685) (0.640) (0.642)
##
## ---------------------------------------------------------------------------------------------
## Observations 526 526 526
## R2 0.165 0.302 0.315
## Adjusted R2 0.163 0.299 0.311
## Residual Std. Error 3.378 (df = 524) 3.092 (df = 523) 3.066 (df = 522)
## F Statistic 103.363*** (df = 1; 524) 113.067*** (df = 2; 523) 79.978*** (df = 3; 522)
## =============================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01stargazer(mod4,mod5,mod_boxcox, type="text", title = "Resultados de regresión 2")#de base, para LaTeX##
## Resultados de regresión 2
## ==============================================================
## Dependent variable:
## -------------------------------
## wage log(wage) wage_boxcox
## (1) (2) (3)
## --------------------------------------------------------------
## educ 0.529*** 0.079*** 0.062***
## (0.047) (0.007) (0.005)
##
## antigüedad 0.154*** 0.019*** 0.015***
## (0.019) (0.003) (0.002)
##
## married 0.683** 0.147*** 0.121***
## (0.275) (0.039) (0.030)
##
## female -1.711*** -0.280*** -0.221***
## (0.266) (0.037) (0.029)
##
## nonwhite -0.065 -0.002 -0.001
## (0.427) (0.060) (0.047)
##
## Constant -1.123* 0.571*** 0.613***
## (0.664) (0.093) (0.073)
##
## --------------------------------------------------------------
## Observations 526 526 526
## R2 0.365 0.400 0.394
## Adjusted R2 0.359 0.394 0.388
## Residual Std. Error (df = 520) 2.956 0.414 0.327
## F Statistic (df = 5; 520) 59.848*** 69.227*** 67.641***
## ==============================================================
## Note: *p<0.1; **p<0.05; ***p<0.016 Modelos lineales generalizados
Al trabajar con datos como los de salarios en el conjunto de datos de Wooldridge, podemos encontrar problemas con las suposiciones fundamentales de los Mínimos Cuadrados Ordinarios (OLS). Aunque los modelos lineales clásicos nos proporcionan resultados, estos pueden ser problemáticos si los errores no están bien distribuidos, es decir, si violan los supuestos de normalidad y homocedasticidad. Esto significa que las inferencias obtenidas pueden no ser adecuadas o estar sesgadas.
Para abordar estas limitaciones, los Modelos Lineales Generalizados (GLM) extienden el marco de los modelos lineales al permitir una mayor flexibilidad en la relación entre las variables independientes y la variable dependiente. Un GLM se compone de tres elementos clave:
- Predictor lineal (\(\eta\)): \(\eta=𝑋𝛽\)
Donde X es la matriz de variables independientes y 𝛽 es el vector de coeficientes.
Función de Enlace (\(𝑔\)): Esta es monótona y diferenciable en todo su dominio, y que transforma el predictor lineal \(g(\mu) = \eta\). Su inversa permite obtener las predicciones de la variable dependiente: \(𝑦̂ = 𝑔⁻¹(𝑋𝛽̂)\), \(\hat{y} = g^{-1}(\eta)\)).
Distribución de Respuesta: La variable dependiente se asume que sigue una distribución de la familia exponencial, denotada por \(f(y | \mu)\) . Algunas distribuciones comunes incluyen la normal, binomial, Poisson y gamma.
Estos componentes proporcionan la flexibilidad necesaria para modelar diferentes tipos de variables dependientes, permitiendo que el modelo capture mejor las características de los datos y aborde problemas como la heterocedasticidad y la no normalidad de los errores, así como cuando el rango de la variable de respuesta está limitado, entre otros.
En R, la función glm() se utiliza para ajustar Modelos Lineales Generalizados. Esta función permite especificar tanto la distribución de la variable dependiente como la función de enlace adecuada.
Sintaxis básica de glm()
glm(formula, family = family_type(link = link_function), data = dataset)formula: Especifica la relación entre las variables dependiente e independientes (similar alm()).family: Define la distribución de la variable dependiente. Puede sergaussian(para regresión lineal),binomial(para regresión logística),poisson(para modelos de conteo),gamma, entre otros.link: Es la función de enlace que conecta la media de la variable dependiente con las variables independientes (por ejemplo,log,identityoinverse).
La función de enlace predeterminada para una familia puede cambiarse especificando un enlace a la función de familia. Si no se informa nada, el modelo correrá en la práctica un modelo lineal simple. Comparen los resultados:
glm(wage ~ educ + antigüedad + married, data = base1)##
## Call: glm(formula = wage ~ educ + antigüedad + married, data = base1)
##
## Coefficients:
## (Intercept) educ antigüedad married
## -2.5258 0.5561 0.1748 0.8924
##
## Degrees of Freedom: 525 Total (i.e. Null); 522 Residual
## Null Deviance: 7160
## Residual Deviance: 4906 AIC: 2677lm(wage ~ educ + antigüedad + married, data = base1)##
## Call:
## lm(formula = wage ~ educ + antigüedad + married, data = base1)
##
## Coefficients:
## (Intercept) educ antigüedad married
## -2.5258 0.5561 0.1748 0.8924Por ejemplo, si la variable de respuesta es no negativa y la varianza es proporcional a la media, se usaría la función de enlace “identity” con la familia “quasipoisson”. Esto se especificaría como:
family = quasipoisson(link = "identity")La decisión sobre qué familia es apropiada no se discute a profundidad en esta secuencia, pero estos pueden ser:
binomial(link = "logit")
gaussian(link = "identity")
Gamma(link = "inverse")
inverse.gaussian(link = "1/mu^2")
poisson(link = "log")
quasi(link = "identity", variance = "constant")
quasibinomial(link = "logit")
quasipoisson(link = "log")- gaussian: Para variables continuas que siguen una distribución normal (equivalente a la regresión lineal clásica).
- binomial: Para variables categóricas binarias (regresión logística).
- poisson: Para datos de conteo (números enteros no negativos).
- Gamma: Para variables continuas y positivas, especialmente cuando la varianza aumenta con la media.
- inverse.gaussian: Para variables continuas positivas con varianza que aumenta rápidamente con la media.
Algunas funciones de enlace comunes son: - identity: Sin transformación (usada en regresión lineal). - log: Transforma la media mediante el logaritmo natural (útil para variables positivas). - logit: Función logística, usada en regresión logística. - inverse: Utiliza la inversa de la media.
Cada familia tiene una función de enlace predeterminada, pero puede
modificarse según las necesidades del análisis. Revisar la página de
ayuda de glm y la documentacion de objetos familiares para
modelos ayudará en gran medida.
Volviendo al ejemplo que traíamos antes, queremos modelar el salario
(wage), que es una variable continua y positiva, y sospechamos que la
varianza aumenta con la media. Podemos utilizar la familia Gamma con una
función de enlace logarítmica, a través de
family = Gamma(link = "log"):
# Ajustando un modelo GLM con distribución gamma y enlace logarítmico
mod_glm <- glm(wage ~ educ + antigüedad + married + female + nonwhite,
family = Gamma(link = "log"), data = base1)
# Resumen del modelo
summary(mod_glm)##
## Call:
## glm(formula = wage ~ educ + antigüedad + married + female +
## nonwhite, family = Gamma(link = "log"), data = base1)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6336436 0.1007043 6.292 6.65e-10 ***
## educ 0.0817722 0.0071560 11.427 < 2e-16 ***
## antigüedad 0.0217897 0.0028409 7.670 8.53e-14 ***
## married 0.1220276 0.0417808 2.921 0.00364 **
## female -0.2856600 0.0404131 -7.068 5.07e-12 ***
## nonwhite -0.0003761 0.0647888 -0.006 0.99537
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2010781)
##
## Null deviance: 158.876 on 525 degrees of freedom
## Residual deviance: 90.711 on 520 degrees of freedom
## AIC: 2305
##
## Number of Fisher Scoring iterations: 6Los coeficientes representan el efecto multiplicativo de las variables independientes sobre el salario. Por pasos:
El intercepto se interpreta como el valor cuando las variables independientes son cero. En este caso, el logaritmo del salario esperado es 0.6336. - Exponenciando el intercepto: \(exp(0.6336) ≈ 1.884\). - Esto significa que, para una persona con cero años de educación y antigüedad, no casada, hombre y blanco, el salario promedio esperado es aproximadamente $1.88 por hora (recordemos, esto es de 1976).
Cada año adicional de educación se asocia con un incremento en el logaritmo del salario de 0.08177. Los coeficientes se han de exponenciar (y son multiplicativos, entendiéndose como cambios porcentuales por unidad adicional).
exp(coef(mod_glm))## (Intercept) educ antigüedad married female nonwhite
## 1.8844643 1.0852086 1.0220288 1.1297853 0.7515181 0.9996239Por ejemplo, por cada año adicional de educación, el salario promedio
aumenta en aproximadamente un 8.53%, manteniendo constantes las demás
variables. En el caso de la variable ‘dummy’, como female, el
exponenciado es menos que uno. Esto implica que las mujeres ganaban, en
promedio, un 24.82% menos que los hombres, manteniendo constantes las
demás variables en el modelo. Otros elementos en las columnas del
estimado son similares a lo que vimos al usar lm().
Sin embargo tenemos información adicional. Hay un parámetro de dispersión. Este es un estimado de la varianza de los residuos en el modelo Gamma. Un valor más pequeño indica menor variabilidad de los datos alrededor del modelo ajustado. La raíz cuadrada de ese parámetro es la desviación estándar estimada. Las desviaciones nulas y residuales son comparaciones entre un modelo vacío (con sólo un intercepto), y el modelo ajustado con todas las variables.
En este caso reporta el criterio de información de Akaike, que es una medida de la calidad del modelo que penaliza por la complejidad (número de parámetros). Al comparar modelos, un AIC más bajo indica un modelo preferido. Como sólo tenemos un modelo, el AIC nos sirve para comparar con futuros modelos alternativos. Finalmente, declara la cantidad de Iteraciones que fueron necesarias para converger y encontrar los estimadores de máxima verosimilitud.
Verifiquemos otra vez con el modelo transformado
# Resumen del modelo
summary(mod5)##
## Call:
## lm(formula = log(wage) ~ educ + antigüedad + married + female +
## nonwhite, data = base1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.89855 -0.27105 -0.03577 0.24487 1.24663
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.570835 0.092937 6.142 1.62e-09 ***
## educ 0.079483 0.006604 12.035 < 2e-16 ***
## antigüedad 0.019443 0.002622 7.416 4.94e-13 ***
## married 0.146703 0.038558 3.805 0.000159 ***
## female -0.280401 0.037296 -7.518 2.45e-13 ***
## nonwhite -0.002419 0.059791 -0.040 0.967738
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4138 on 520 degrees of freedom
## Multiple R-squared: 0.3996, Adjusted R-squared: 0.3939
## F-statistic: 69.23 on 5 and 520 DF, p-value: < 2.2e-16Notamos que la transformación aproxima los coeficientes a los
conseguidos con el glm().
library(arm)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loading required package: lme4##
## arm (Version 1.14-4, built: 2024-4-1)## Working directory is /Users/rashid/Library/CloudStorage/OneDrive-UniversityofPuertoRico/Universidad de Puerto Rico/Servicio/Secuencia-de-Talleres-RStudio-FCS-UPRRP-Sept-2025/Taller 3##
## Attaching package: 'arm'## The following object is masked from 'package:car':
##
## logit## The following object is masked from 'package:jtools':
##
## standardize# Gráfico de coeficientes para el modelo GLM
coefplot(mod_glm, col.pts = "red", cex.pts = 1.5, main = "Comparación de Coeficientes")
# Añadimos los coeficientes del modelo OLS al mismo gráfico
coefplot(mod5, add = TRUE, col.pts = "blue", cex.pts = 1.5)
legend("topright", legend = c("GLM (Gamma con enlace log)", "OLS (log(wage))"),
col = c("red", "blue"), pch = 16)
6.1 Análisis de la bondad del modelo
# Generamos los gráficos diagnósticos
par(mfrow = c(2, 2))
plot(mod_glm)
Esto se ve posiblemente mejor que la variante de OLS, aunque tampoco parece perfecto.
7 Modelos jerárquicos
Los modelos jerárquicos, también conocidos como modelos lineales mixtos o modelos multinivel, son una extensión de los modelos lineales que permiten analizar datos en los que las observaciones están agrupadas o anidadas en diferentes niveles. Estos modelos son especialmente útiles cuando se espera que exista correlación entre las observaciones dentro de los mismos grupos.
En este análisis, utilizaremos el paquete lme4 en R, que
es ampliamente utilizado para ajustar modelos lineales mixtos. A
continuación, te guiaré a través de los pasos para ajustar y comprender
estos modelos, explicando cada componente y resultado de manera clara.
La sintaxis de lme4 se basa en la sintaxis de los modelos
lineales que ya conocemos de lm().
La función lmer() de lme4 añade la
especificación de la variable de grupo/sujeto y de la estructura de
efectos aleatorios que se van a estimar. En paréntesis adicionales (ver
abajo), el término a la izquierda de | especifica los
efectos aleatorios que se van a estimar. El término a la derecha de
| representa la(s) variable(s) que definen la estructura de
agrupación (o anidamiento) de los datos.
Un 1 en la parte izquierda del paréntesis significa que
se debe estimar un componente de varianza de intercepto aleatorio. Si
también se coloca la variable predictora a la izquierda de
|, esto indica que se deben incluir pendientes aleatorias.
La forma básica de estos será:
lmer(data = datos, VarDependiente ~ VarIndependiente + (1 | VarGrupal))Corresponde a un modelo que puede describirse de la siguiente manera: “La variable dependiente es predicha por la variable independiente. Al mismo tiempo, la varianza de los residuos de nivel 2 del intercepto es un parámetro del modelo”.
Comencemos con modelos HLM que incluyen solo variables predictoras de nivel 1 (individuos anidados en grupos). Luego, en un segundo paso, añadiremos predictores de nivel 2. Finalmente, analizaremos los modelos HLM de medidas repetidas, donde el nivel más bajo (nivel 1) corresponde a observaciones repetidas dentro de los individuos, y esas observaciones están anidadas en los participantes individuales (nivel 2).
Nos limitaremos aquí a modelos de 2 niveles. Los principios de los modelos HLM pueden ilustrarse de manera bastante parsimoniosa de esta forma, y expandir los modelos a más de dos niveles de análisis es bastante sencillo.
# Cargar las librerías necesarias
#library(tidyverse)
library(lme4)
library(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## step# Cargar los datos desde una URL
df <- read_csv("https://raw.githubusercontent.com/methodenlehre/data/master/salary-data.csv")## Rows: 600 Columns: 4## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): firma, sector
## dbl (2): experience, salary
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# Convertir las variables 'firma' y 'sector' en factores
df <- df |>
mutate(firma = as.factor(firma),
sector = as.factor(sector))
# Ver las últimas filas del conjunto de datos
tail(df)## # A tibble: 6 × 4
## firma experience salary sector
## <fct> <dbl> <dbl> <fct>
## 1 Firma 20 3.58 6838. Privat
## 2 Firma 20 3.18 7604. Privat
## 3 Firma 20 3.39 5714. Privat
## 4 Firma 20 7.12 10089. Privat
## 5 Firma 20 2.98 6940. Privat
## 6 Firma 20 6.45 9330. Privat7.1 Modelo nulo
En este caso saco unos datos de salarios de compañías ficticias en
Suiza. Con mutate() convertí las variables firma (empresa)
y sector en factores, ya que representan categorías. Antes de añadir
predictores, ajustamos un modelo nulo que solo incluye el intercepto y
un término aleatorio para capturar la variabilidad entre las empresas
(firma). Este modelo nos permite estimar la correlación intraclase (ICC)
y entender qué proporción de la variabilidad total del salario se debe a
diferencias entre empresas.
library(Matrix)
library(lme4)
# Ajustar el modelo nulo con intercepto aleatorio por empresa
modelo_nulo <- lmer(salary ~ 1 + (1 | firma), data = df, REML = TRUE)
# Obtener las predicciones del modelo nulo
df$predicciones_nulo <- predict(modelo_nulo)
# Resumen del modelo
summary(modelo_nulo)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: salary ~ 1 + (1 | firma)
## Data: df
##
## REML criterion at convergence: 10433.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9816 -0.6506 -0.0494 0.5779 4.2131
##
## Random effects:
## Groups Name Variance Std.Dev.
## firma (Intercept) 851249 922.6
## Residual 1954745 1398.1
## Number of obs: 600, groups: firma, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8737.6 214.1 19.0 40.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Explicación:
- lmer(): Ajusta un modelo lineal mixto.
salary ~ 1: Indica que solo se ajusta el intercepto fijo (media general del salario).(1 | firma): Especifica un intercepto aleatorio para cada empresa (firma), capturando la variabilidad entre empresas.REML = TRUE: Utiliza el método de Máxima Verosimilitud Restringida para estimar los parámetros, lo cual es apropiado cuando se comparan modelos con diferentes efectos fijos.predict(): Genera predicciones del modelo para cada observación.
La varianza del intercepto aleatorio (firma) representa la variabilidad del salario promedio entre empresas. Por otro lado la varianza residual captura la variabilidad del salario dentro de las empresas. Calcularemos la correlación intraclase como la proporción de la varianza total que se debe a diferencias entre empresas. Se define como: \(\rho = \frac{\sigma^2_{\text{Nivel-2}}}{\sigma^2_{\text{Nivel-2}} + \sigma^2_{\text{Nivel-1}}}\)
# Extraer las varianzas del modelo
var_intercepto <- as.numeric(VarCorr(modelo_nulo)$firma[1])
var_residual <- attr(VarCorr(modelo_nulo), "sc")^2
# Calcular la ICC
ICC <- var_intercepto / (var_intercepto + var_residual)
ICC## [1] 0.303368Correlación intraclase: \(\hat{\rho} = \frac{\hat{\sigma}^2_{υ_0}}{\hat{\sigma}^2_{υ_0} + \hat{\sigma}^2_{ϵ}} = \frac{851249}{851249 + 1954745} = 0.3034\)
Es decir, 30% de la variabilidad total del salario se debe a diferencias entre empresas. Existe una correlación notable entre los salarios de los empleados dentro de la misma empresa, lo que justifica el uso de un modelo multinivel.
Con la función ranef() podemos ver los efectos aleatorios (residuos de segundo nivel del intercepto):
ranef(modelo_nulo)## $firma
## (Intercept)
## Firma 01 789.65416
## Firma 02 1105.01807
## Firma 03 1923.02618
## Firma 04 -1136.50080
## Firma 05 953.54595
## Firma 06 -958.93185
## Firma 07 -10.14627
## Firma 08 -254.84048
## Firma 09 -651.31802
## Firma 10 768.48590
## Firma 11 -506.26620
## Firma 12 940.40709
## Firma 13 -742.84383
## Firma 14 -975.45482
## Firma 15 -1161.36931
## Firma 16 -97.68008
## Firma 17 661.96052
## Firma 18 -168.11195
## Firma 19 351.23926
## Firma 20 -829.87351
##
## with conditional variances for "firma"Para determinar si la varianza entre empresas es estadísticamente significativa, comparamos el modelo nulo con un modelo sin efectos aleatorios (modelo lineal simple).
# Utilizar la función ranova() para comparar modelos
anova_aleatorio <- ranova(modelo_nulo)
anova_aleatorio## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## salary ~ (1 | firma)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 3 -5216.7 10440
## (1 | firma) 2 -5295.5 10595 157.44 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Este comparó el modelo actual con un modelo reducido que elimina el efecto aleatorio. El valor p indica que incluir el intercepto aleatorio para firma mejora significativamente el modelo. Aproximadamente, el 30% de la varianza total del salario se puede atribuir a diferencias entre empresas.
7.2 Modelo con predictor de primer nivel
Ahora, añadimos un predictor de nivel 1 (experience), que es una variable individual, al modelo. Mantenemos el intercepto aleatorio para capturar la variabilidad entre empresas.
# Ajustar el modelo con 'experience' como predictor fijo y intercepto aleatorio por empresa
modelo_intercepto_aleatorio <- lmer(salary ~ experience + (1 | firma), data = df, REML = TRUE)
# Obtener las predicciones del modelo
df$predicciones_intercepto_aleatorio <- predict(modelo_intercepto_aleatorio)
# Resumen del modelo
summary(modelo_intercepto_aleatorio)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: salary ~ experience + (1 | firma)
## Data: df
##
## REML criterion at convergence: 10127.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8109 -0.6884 0.0005 0.5980 3.8833
##
## Random effects:
## Groups Name Variance Std.Dev.
## firma (Intercept) 614367 783.8
## Residual 1184502 1088.3
## Number of obs: 600, groups: firma, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5964.18 229.41 48.00 26.00 <2e-16 ***
## experience 534.34 27.21 589.48 19.64 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## experience -0.615Los resultados del modelo mixto muestran que el efecto fijo de la experiencia es significativamente positivo, con un coeficiente estimado de \(\hat{\gamma}_{10} = 534.34\) y un valor \(p\) menor a 0.001. Esto indica que, en promedio, el salario aumenta en aproximadamente 534 unidades (por ejemplo, francos suizos) por cada año adicional de experiencia, controlando por las diferencias entre empresas. El modelo incluye un intercepto aleatorio a nivel de empresa (firma), lo que nos permite ajustar diferencias entre las empresas en los salarios base.
Se le puede añadir un nivel adicional al modelo para tanto acomodar por un intercepto así como pendiente aleatoria. Al aplicar este modelo la varianza (aleatoria) de la pendiente de la variable independiente de primer nivel se puede estimar, dándonos una idea de las diferencias en el efecto (pendiente) de experiencia en salario entre las unidades del segundo nivel (firmas). Si comparamos los efectos entre este y el modelo anterior, vemos la disminución marginal de la varianza del intercepto aleatorio en comparación con el modelo nulo, lo que sugiere que parte de la variabilidad entre empresas se explica por la experiencia de los empleados.
fixef(modelo_intercepto_aleatorio)## (Intercept) experience
## 5964.1757 534.3446ranef(modelo_intercepto_aleatorio)## $firma
## (Intercept)
## Firma 01 204.30371
## Firma 02 646.14732
## Firma 03 1492.00151
## Firma 04 -910.78990
## Firma 05 389.16512
## Firma 06 -924.63977
## Firma 07 577.66959
## Firma 08 -516.51767
## Firma 09 -638.24646
## Firma 10 768.48113
## Firma 11 -619.55111
## Firma 12 1091.33530
## Firma 13 -773.67207
## Firma 14 -738.17926
## Firma 15 -652.94087
## Firma 16 57.33923
## Firma 17 458.05487
## Firma 18 -89.38416
## Firma 19 944.22822
## Firma 20 -764.80474
##
## with conditional variances for "firma"7.3 Modelo con pendientes aleatorias
Para investigar si el efecto de la experiencia en el salario varía entre empresas, ajustamos un modelo que incluye tanto intercepto como pendiente aleatorios para experience.
# Ajustar el modelo con intercepto y pendiente aleatorios por empresa
modelo_coeficientes_aleatorios <- lmer(salary ~ experience + (experience | firma), data = df, REML = TRUE)
# Obtener las predicciones del modelo
df$predicciones_coeficientes_aleatorios <- predict(modelo_coeficientes_aleatorios)
# Resumen del modelo
summary(modelo_coeficientes_aleatorios)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: salary ~ experience + (experience | firma)
## Data: df
##
## REML criterion at convergence: 10117.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8307 -0.6804 0.0037 0.5999 3.3608
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## firma (Intercept) 722588 850.1
## experience 18392 135.6 -0.51
## Residual 1136296 1066.0
## Number of obs: 600, groups: firma, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5933.72 240.38 18.89 24.68 7.77e-16 ***
## experience 530.85 40.59 18.95 13.08 6.20e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## experience -0.690En este caso (experience | firma) especifica que tanto
el intercepto como la pendiente de experience varían aleatoriamente
entre empresas.
Los resultados muestran que el efecto fijo de la experiencia sigue siendo significativamente positivo, con un coeficiente estimado de \(\hat{\gamma}_{10} = 530.85\) y un valor \(p\) menor a 0.001. Esto indica que, en promedio, el salario aumenta en aproximadamente 531 unidades (por ejemplo, CHF para estos datos ficticios) por cada año adicional de experiencia. Una covarianza negativa sugiere que las empresas con salarios base más altos tienden a tener incrementos más pequeños por año de experiencia, y viceversa.
Para evaluar si el modelo con pendientes aleatorias proporciona un mejor ajuste que el modelo con solo intercepto aleatorio, realizamos una comparación de modelos.
# Ajustar el modelo con intercepto aleatorio (modelo reducido)
modelo_intercepto_aleatorio <- lmer(salary ~ experience + (1 | firma), data = df, REML = FALSE)
# Ajustar el modelo con intercepto y pendiente aleatorios (modelo completo)
modelo_coeficientes_aleatorios <- lmer(salary ~ experience + (experience | firma), data = df, REML = FALSE)
# Comparar los modelos utilizando ANOVA
anova(modelo_intercepto_aleatorio, modelo_coeficientes_aleatorios)## Data: df
## Models:
## modelo_intercepto_aleatorio: salary ~ experience + (1 | firma)
## modelo_coeficientes_aleatorios: salary ~ experience + (experience | firma)
## npar AIC BIC logLik -2*log(L) Chisq Df
## modelo_intercepto_aleatorio 4 10156 10173 -5073.9 10148
## modelo_coeficientes_aleatorios 6 10150 10177 -5069.2 10138 9.3127 2
## Pr(>Chisq)
## modelo_intercepto_aleatorio
## modelo_coeficientes_aleatorios 0.009501 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El p-valor pequeño indica que el modelo con pendientes aleatorias es significativamente mejor que el modelo con solo intercepto aleatorio. Existe variación significativa en el efecto de la experiencia en el salario entre empresas.
Para ayudar a comprender mejor los resultados, es útil visualizar cómo varía el efecto de la experiencia en el salario entre empresas.
#library(ggplot2)
# Extraer los coeficientes aleatorios
coeficientes_aleatorios <- coef(modelo_coeficientes_aleatorios)$firma
# Renombrar las columnas
colnames(coeficientes_aleatorios) <- c("Intercepto", "Pendiente")
# Añadir el identificador de la empresa
coeficientes_aleatorios$firma <- rownames(coeficientes_aleatorios)
# Graficar las pendientes de experiencia por empresa
ggplot(coeficientes_aleatorios, aes(x = reorder(firma, Pendiente), y = Pendiente)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(title = "Pendientes de Experiencia por Empresa",
x = "Empresa",
y = "Efecto de la Experiencia en el Salario") +
theme_minimal()
# Crear un conjunto de datos con las predicciones del modelo
df$predicciones <- predict(modelo_coeficientes_aleatorios)
# Graficar los datos y las líneas de regresión por empresa
ggplot(df, aes(x = experience, y = salary, color = firma, group = firma)) +
geom_point(alpha = 0.5) +
geom_line(aes(y = predicciones)) +
labs(title = "Relación entre experiencia y salario por empresa",
x = "Experiencia",
y = "Salario") +
theme_minimal()
8 Modelos longitudinales
Los modelos longitudinales son una herramienta estadística esencial para analizar datos en los que se realizan observaciones repetidas de los mismos individuos a lo largo del tiempo. Este tipo de análisis es fundamental en campos como la medicina, psicología, economía, política, sociología y otras disciplinas donde es importante entender cómo cambian las mediciones en los individuos o unidades a lo largo del tiempo y qué factores influyen en esos cambios.
Los datos longitudinales son un caso especial de datos jerárquicos o multinivel, donde las observaciones (mediciones) están anidadas dentro de individuos. En este contexto, las observaciones repetidas de un mismo individuo tienden a ser más similares entre sí que las observaciones de diferentes individuos. Esta correlación dentro de los sujetos es crucial y debe ser tenida en cuenta para obtener inferencias estadísticas válidas. Una vista general sobre paquetes útiles en R para analizar este tipo de datos correlacionados se puede encontrar en la CRAN Task View dedicada.
En este taller, exploraremos cómo ajustar y analizar modelos longitudinales utilizando R, centrándonos en un conjunto de datos que mide el crecimiento dental en niños y niñas a diferentes edades. Trabajaremos con un conjunto de datos que contiene medidas longitudinales del crecimiento dental en 27 niños (16 niños y 11 niñas). Las medidas corresponden a la distancia pituitaria-pterigomaxilar (una medida de crecimiento dental) y se tomaron a las edades de 8, 10, 12 y 14 años. Nuestro objetivo es describir cómo cambia esta distancia con la edad y comparar el patrón de crecimiento entre niños y niñas.
Primero, cargamos los datos y observamos su estructura:
load(url("http://alecri.github.io/downloads/data/dental.RData"))
head(dental)## # A tibble: 6 × 6
## id sex y8 y10 y12 y14
## <dbl> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 1 Girl 21 20 21.5 23
## 2 2 Girl 21 21.5 24 25.5
## 3 3 Girl 20.5 24 24.5 26
## 4 4 Girl 23.5 24.5 25 26.5
## 5 5 Girl 21.5 23 22.5 23.5
## 6 6 Girl 20 21 21 22.5Los datos se presentan en formato ancho, donde las mediciones repetidas se encuentran en columnas separadas para cada edad. Este formato no es ideal para el análisis longitudinal, por lo que convertiremos los datos a un formato largo usando la función pivot_longer():
library(labelled) # etiquetado de datos
library(rstatix) # estadísticas descriptivas##
## Attaching package: 'rstatix'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:stats':
##
## filterlibrary(ggpubr) # estadísticas descriptivas y gráficos convenientes
library(GGally) # gráficos avanzados
library(car) # útil para ANOVA/pruebas de Wald
library(Epi) # fácil obtención de intervalos de confianza para coeficientes/predicciones del modelo##
## Attaching package: 'Epi'## The following object is masked from 'package:lme4':
##
## factorize#library(lme4) # modelos lineales de efectos mixtos
#library(lmerTest) # pruebas para modelos lineales de efectos mixtos
library(emmeans) # medias marginales## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'##
## Attaching package: 'emmeans'## The following object is masked from 'package:GGally':
##
## pigslibrary(multcomp) # intervalos de confianza para combinaciones lineales de coeficientes del modelo## Loading required package: mvtnorm##
## Attaching package: 'mvtnorm'## The following object is masked from 'package:arm':
##
## standardize## The following object is masked from 'package:jtools':
##
## standardize## Loading required package: survival## Loading required package: TH.data##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyserlibrary(geepack) # ecuaciones de estimación generalizadas
library(ggeffects) # efectos marginales, predicciones ajustadas##
## Attaching package: 'ggeffects'## The following object is masked from 'package:wooldridge':
##
## fishlibrary(gt) # tablas bonitas
dental_long <- pivot_longer(dental, cols = starts_with("y"),
names_to = "measurement", values_to = "distance") |>
mutate(
age = parse_number(measurement),
measurement = fct_inorder(paste("Medida a los", age))
) |>
set_variable_labels(
age = "Edad del niño/a al momento de la medición",
measurement = "Etiqueta de medición temporal",
distance = "Medición de distancia"
)
head(dental_long)## # A tibble: 6 × 5
## id sex measurement distance age
## <dbl> <fct> <fct> <dbl> <dbl>
## 1 1 Girl Medida a los 8 21 8
## 2 1 Girl Medida a los 10 20 10
## 3 1 Girl Medida a los 12 21.5 12
## 4 1 Girl Medida a los 14 23 14
## 5 2 Girl Medida a los 8 21 8
## 6 2 Girl Medida a los 10 21.5 10Ahora, cada fila representa una medición de un individuo en una edad específica, lo que facilita el análisis longitudinal. Antes de ajustar modelos, exploramos los datos de forma descriptiva para entender las tendencias generales.
group_by(dental_long, age) |>
get_summary_stats(distance)## # A tibble: 4 × 14
## age variable n min max median q1 q3 iqr mad mean sd
## <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 8 distance 27 16.5 27.5 22 21 23.2 2.25 1.48 22.2 2.43
## 2 10 distance 27 19 28 23 21.5 24.5 3 2.22 23.2 2.16
## 3 12 distance 27 19 31 24 23 26 3 2.22 24.6 2.82
## 4 14 distance 27 19.5 31.5 26 25 27.8 2.75 2.22 26.1 2.77
## # ℹ 2 more variables: se <dbl>, ci <dbl>Observamos que la distancia media aumenta con la edad, lo que sugiere un crecimiento dental a medida que los niños y niñas crecen. Creamos un diagrama de caja para visualizar la distribución de la distancia a cada edad:
ggplot(dental_long, aes(measurement, distance, fill = measurement)) +
geom_boxplot() +
geom_jitter(width = 0.2) +
guides(fill = "none") +
labs(x = "Edad (años)", y = "Crecimiento dental, mm")
El gráfico muestra que la mediana y los valores de la distancia aumentan con la edad. La dispersión también parece aumentar ligeramente, lo que indica mayor variabilidad en edades superiores. Exploremos si hay diferencias en el crecimiento dental entre niños y niñas.
group_by(dental_long, sex, measurement) |>
get_summary_stats(distance, show = c("mean", "sd"))## # A tibble: 8 × 6
## sex measurement variable n mean sd
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Girl Medida a los 8 distance 11 21.2 2.12
## 2 Girl Medida a los 10 distance 11 22.2 1.90
## 3 Girl Medida a los 12 distance 11 23.1 2.37
## 4 Girl Medida a los 14 distance 11 24.1 2.44
## 5 Boy Medida a los 8 distance 16 22.9 2.45
## 6 Boy Medida a los 10 distance 16 23.8 2.14
## 7 Boy Medida a los 12 distance 16 25.7 2.65
## 8 Boy Medida a los 14 distance 16 27.5 2.08Y gráficamente
ggplot(dental_long, aes(sex, distance, fill = measurement)) +
geom_boxplot() +
labs(x = "", y = "Crecimiento dental, mm", fill = "")
Visualmente, parece que los niños tienen, en promedio, una distancia ligeramente mayor que las niñas en cada edad.
Evaluemos las correlaciones entre medidas de interés:
ggpairs(select(dental, starts_with("y")), lower = list(continuous = "smooth"))
y separemos nuevamente por grupo (sexo)
ggpairs(dental, mapping = aes(colour = sex,alpha = 0.5), columns = 3:6,
lower = list(continuous = "smooth"))
Digamos que queremos entender si el efecto varía con el tiempo. Podríamos evaluar gráficamente esto:
group_by(dental_long, sex, age) |>
summarise(mean = list(mean_ci(distance)), .groups = "drop") |>
unnest_wider(mean) |>
mutate(agex = age - .05 + .05*(sex == "Boy")) |>
ggplot(aes(agex, y, col = sex, shape = sex)) +
geom_point() +
geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.2) +
geom_line() +
labs(x = "Edad, en años", y = "Crecimiento dental medio, mm", shape = "Sex", col = "Sex")
ggplot(dental_long, aes(x = age, y = distance, color = sex)) +
geom_point() +
geom_line(aes(group = id)) +
labs(x = "Edad (años)", y = "Crecimiento dental (mm)", color = "Sexo") +
theme_minimal()
Para analizar los datos longitudinales, utilizamos modelos mixtos lineales, que permiten modelar tanto los efectos fijos (comunes a todos los individuos) como los efectos aleatorios (específicos de cada individuo).
8.1 Modelo inicial
Ajustamos un modelo que considera la edad como factor categórico (variables indicadoras) y un intercepto aleatorio para capturar la variabilidad entre individuos.
# Ajustar el modelo mixto
modelo_edad <- lmer(distance ~ factor(age) + (1 | id), data = dental_long)
# Resumen del modelo
summary(modelo_edad)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: distance ~ factor(age) + (1 | id)
## Data: dental_long
##
## REML criterion at convergence: 443.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7376 -0.5248 0.0153 0.4027 3.7212
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 4.465 2.113
## Residual 2.078 1.442
## Number of obs: 108, groups: id, 27
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 22.1852 0.4923 43.3911 45.066 < 2e-16 ***
## factor(age)10 0.9815 0.3924 78.0000 2.501 0.0145 *
## factor(age)12 2.4630 0.3924 78.0000 6.277 1.80e-08 ***
## factor(age)14 3.9074 0.3924 78.0000 9.958 1.52e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) fc()10 fc()12
## factor(g)10 -0.399
## factor(g)12 -0.399 0.500
## factor(g)14 -0.399 0.500 0.500En este caso el intercepto captura la distancia media dental a los 8 años (valor de referencia). Esta es 22.18 mm. La diferencia entre las respuestas medias de los niños de 10 y 8 años es de 0.98 mm. La diferencia entre las respuestas medias de los niños de 12 y 8 años es de 2.46 mm. La diferencia entre las respuestas medias de los niños de 14 y 8 años es de 3.91 mm.
Esto indica que la distancia pituitaria-pterigomaxilar aumenta significativamente con la edad, lo que es consistente con el crecimiento observado durante este período.
Usamos un análisis de varianza para evaluar si la edad tiene un efecto significativo en el crecimiento dental.
library(car)
Anova(modelo_edad)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: distance
## Chisq Df Pr(>Chisq)
## factor(age) 114.12 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El efecto de la edad es altamente significativo (p < 0.001), indicando que la edad influye en el crecimiento dental.
Calculamos las medias marginales estimadas para cada edad y sus intervalos de confianza.
library(emmeans)
# Calcular las medias marginales
medias_edad <- emmeans(modelo_edad, ~ age)
# Mostrar los resultados
summary(medias_edad)## age emmean SE df lower.CL upper.CL
## 8 22.2 0.492 43.4 21.2 23.2
## 10 23.2 0.492 43.4 22.2 24.2
## 12 24.6 0.492 43.4 23.7 25.6
## 14 26.1 0.492 43.4 25.1 27.1
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95Las medias estimadas confirman que la distancia media aumenta con la edad, aumentando la tasa de crecimiento en cada momento de medición. Podemos presentar la trayectoria estimada del modelo ajustado utilizando la función predict:
# Generar predicciones del modelo ajustado
predicciones <- predict(modelo_edad)
# Visualizar la trayectoria estimada
ggplot(dental_long, aes(x = age, y = distance, color = sex)) +
geom_point() +
geom_line(aes(y = predicciones, group = id), linetype = "dashed") +
labs(x = "Edad (años)", y = "Distancia pituitaria-pterigomaxilar (mm)", color = "Sexo") +
theme_minimal()
Este gráfico muestra tanto los datos reales de distancia en función de la edad como las predicciones del modelo para cada niño, con líneas discontinuas que representan las trayectorias estimadas a lo largo del tiempo.
8.2 Modelo con variable de sexo
# Ajustar el modelo con sexo y su interacción con la edad
modelo_sexo <- lmer(distance ~ factor(age) * sex + (1 | id), data = dental_long)
# Resumen del modelo
summary(modelo_sexo)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: distance ~ factor(age) * sex + (1 | id)
## Data: dental_long
##
## REML criterion at convergence: 423.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7358 -0.4761 0.0469 0.4691 3.6613
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 3.285 1.813
## Residual 1.975 1.405
## Number of obs: 108, groups: id, 27
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 21.1818 0.6915 46.0791 30.630 < 2e-16 ***
## factor(age)10 1.0455 0.5992 75.0000 1.745 0.0851 .
## factor(age)12 1.9091 0.5992 75.0000 3.186 0.0021 **
## factor(age)14 2.9091 0.5992 75.0000 4.855 6.41e-06 ***
## sexBoy 1.6932 0.8983 46.0791 1.885 0.0658 .
## factor(age)10:sexBoy -0.1080 0.7784 75.0000 -0.139 0.8901
## factor(age)12:sexBoy 0.9347 0.7784 75.0000 1.201 0.2337
## factor(age)14:sexBoy 1.6847 0.7784 75.0000 2.164 0.0336 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) fc()10 fc()12 fc()14 sexBoy f()10: f()12:
## factor(g)10 -0.433
## factor(g)12 -0.433 0.500
## factor(g)14 -0.433 0.500 0.500
## sexBoy -0.770 0.334 0.334 0.334
## fctr(g)10:B 0.334 -0.770 -0.385 -0.385 -0.433
## fctr(g)12:B 0.334 -0.385 -0.770 -0.385 -0.433 0.500
## fctr(g)14:B 0.334 -0.385 -0.385 -0.770 -0.433 0.500 0.500Interpretando salidas:
La varianza del Intercepto Aleatorio (id) es 3.285 con una desviación estándar de 1.813. Existe variabilidad significativa entre los individuos en sus medidas iniciales de distancia dental.
La varianza residual es 1.975 con una desviación estándar de 1.405. Esta captura variabilidad no explicadas por el modelo.
El intercepto captura el ser una niña de 8 años como el punto de referencia: 21.1818mm.
Coeficientes de factor(g)AA: El efecto de edad está dándonos el cambio en la distancia dental en comparación con la edad de referencia (8 años) para las niñas. Por ejemplo, las niñas de 10 años tienen, en promedio, una distancia dental 1.05 mm mayor que las niñas de 8 años, mientras que las de 14 años tienen, en promedio, una distancia dental 2.91 mm mayor que las niñas de 8 años.
Coeficiente
sexBoy: Los niños de 8 años tienen, en promedio, una distancia dental 1.69 mm mayor que las niñas de 8 años. Esta es una variable dummy.Los coeficientes de interacción indican cuánto cambia la diferencia entre niños y niñas en cada edad en comparación con la edad de referencia (8 años). Por ejemplo, a los 14 años, los niños tienen una distancia dental adicional de 1.68 mm en comparación con las niñas, más allá de la diferencia observada a los 8 años. Esa edad es la única que reporta distancias significativas a la referencia; al entrar más en la adolescencia las diferencias se tornan notables.
Anova(modelo_sexo)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: distance
## Chisq Df Pr(>Chisq)
## factor(age) 120.0950 3 < 2.2e-16 ***
## sex 9.2921 1 0.002301 **
## factor(age):sex 7.0847 3 0.069247 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El efecto de la edad es altamente significativo. Esto indica que existen diferencias significativas en la distancia dental media entre las diferentes edades (8, 10, 12 y 14 años), independientemente del sexo. También el sexo del paciente parece ser significativo. La interacción entre edad y sexo es marginalmente significativa (p < 0.10, pero > 0.05).
9 Series de tiempo: ejemplo de finanzas
Las series de tiempo son un conjunto de observaciones registradas en momentos sucesivos en el tiempo, ordenadas cronológicamente. El análisis de series de tiempo es fundamental en diversos campos como la economía, finanzas, meteorología, y ciencias sociales, ya que permite comprender y predecir comportamientos futuros basados en datos históricos.
En esta ñapita del taller, me enfoco en el análisis de series de tiempo financieras, específicamente en el estudio de los precios de las acciones. El objetivo es proporcionar una guía clara y práctica sobre modelos de series de tiempo utilizando R.
Las series de tiempo financieras, como los precios de las acciones o los tipos de cambio, son intrínsecamente volátiles y están influenciadas por múltiples factores económicos y políticos. El análisis de estas series permite a los inversores y analistas comprender las tendencias del mercado, evaluar riesgos, anticiparse a estas y tomar decisiones informadas.
Para este análisis, utilizaremos datos históricos del precio de cierre ajustado de las acciones de Apple Inc. (AAPL). Usaremos el paquete quantmod para descargar los datos directamente desde Yahoo Finance.
# Instalar paquetes si no están instalados
#install.packages("quantmod")
#install.packages("forecast")
#install.packages("tseries")
#install.packages("ggplot2")
#install.packages("quantmod", repos = "https://cloud.r-project.org/")
# Cargar los paquetes
library(quantmod)## Loading required package: xts##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## #################################################################################
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, last## Loading required package: TTR##
## Attaching package: 'TTR'## The following object is masked from 'package:Epi':
##
## ROC## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(forecast)##
## Attaching package: 'forecast'## The following object is masked from 'package:ggpubr':
##
## gghistogramlibrary(tseries)
library(ggplot2)Descargamos los datos de los últimos cinco años.
# Establecer el rango de fechas
fecha_inicio <- as.Date("2019-01-01")
fecha_fin <- as.Date(Sys.Date())
# Descargar los datos de AAPL
getSymbols("AAPL", src = "yahoo", from = fecha_inicio, to = fecha_fin)## [1] "AAPL"# Ver las primeras filas
head(AAPL)## AAPL.Open AAPL.High AAPL.Low AAPL.Close AAPL.Volume AAPL.Adjusted
## 2019-01-02 38.7225 39.7125 38.5575 39.4800 148158800 37.57521
## 2019-01-03 35.9950 36.4300 35.5000 35.5475 365248800 33.83243
## 2019-01-04 36.1325 37.1375 35.9500 37.0650 234428400 35.27673
## 2019-01-07 37.1750 37.2075 36.4750 36.9825 219111200 35.19820
## 2019-01-08 37.3900 37.9550 37.1300 37.6875 164101200 35.86919
## 2019-01-09 37.8225 38.6325 37.4075 38.3275 180396400 36.47833Extraemos la columna de precios de cierre ajustados y creamos un objeto de serie de tiempo.
# Extraer el precio de cierre ajustado
precio_cierre <- AAPL[, "AAPL.Adjusted"]
# Convertir a serie de tiempo
serie_tiempo <- ts(precio_cierre, frequency = 252) # 252 días hábiles en un añoAntes de ajustar cualquier modelo, es importante entender las características de la serie de tiempo.
# Graficar la serie de tiempo
autoplot(serie_tiempo) +
labs(title = "Precio de Cierre Ajustado de AAPL",
x = "Tiempo",
y = "Precio ($)") +
theme_minimal()
El gráfico muestra la evolución del precio de AAPL desde 2019 hasta la fecha actual. Observamos tendencias ascendentes y periodos de volatilidad, especialmente durante eventos económicos significativos como la pandemia de COVID-19 y la subsiguiente alza en valuaciones financieras.
Descomponemos la serie para analizar sus componentes: tendencia, estacionalidad y residuales.
# Verificar la estructura de 'serie_tiempo'
str(serie_tiempo)## Time-Series [1:1703, 1] from 1 to 7.75: 37.6 33.8 35.3 35.2 35.9 ...
## - attr(*, "index")= num [1:1703] 1.55e+09 1.55e+09 1.55e+09 1.55e+09 1.55e+09 ...
## ..- attr(*, "tzone")= chr "UTC"
## ..- attr(*, "tclass")= chr "Date"
## - attr(*, "src")= chr "yahoo"
## - attr(*, "updated")= POSIXct[1:1], format: "2025-10-10 12:49:04"
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr "AAPL.Adjusted"# Convertir a vector numérico
precio_cierre_vector <- as.numeric(precio_cierre)
# Crear la serie de tiempo
serie_tiempo <- ts(precio_cierre_vector, frequency = 252)
# Descomposición utilizando STL (Seasonal and Trend decomposition using Loess)
descomposición <- stl(serie_tiempo, s.window = "periodic")
# Graficar la descomposición
autoplot(descomposición) +
labs(title = "Descomposición de la Serie de Tiempo de AAPL")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the forecast package.
## Please report the issue at <https://github.com/robjhyndman/forecast/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
- Tendencia: Muestra el movimiento a largo plazo del precio.
- Estacionalidad: En series financieras diarias, la estacionalidad puede no ser pronunciada, pero pueden existir patrones semanales o mensuales.
- Residuales: Parte de la serie no explicada por la tendencia ni la estacionalidad.
La estacionariedad es una propiedad clave en el análisis de series de tiempo. Una serie estacionaria tiene estadísticas (media, varianza) constantes en el tiempo. Realizamos la prueba Dickey-Fuller Aumentada (ADF) para verificar si la serie es estacionaria.
# Prueba ADF
adf.test(serie_tiempo, alternative = "stationary")##
## Augmented Dickey-Fuller Test
##
## data: serie_tiempo
## Dickey-Fuller = -2.9984, Lag order = 11, p-value = 0.1557
## alternative hypothesis: stationaryEl p-valor es alto (mayor que 0.05), lo que indica que no podemos rechazar la hipótesis nula de que la serie tiene una raíz unitaria (no estacionaria). Concluimos que la serie no es estacionaria.
Para lograr estacionariedad, transformamos los precios en rendimientos logarítmicos.
# Calcular los rendimientos logarítmicos
rendimientos <- diff(log(serie_tiempo))
# Graficar los rendimientos
autoplot(rendimientos) +
labs(title = "Rendimientos Logarítmicos Diarios de AAPL",
x = "Tiempo",
y = "Rendimiento") +
theme_minimal()
Verificamos nuevamente la serie transformada con ADF.
# Prueba ADF en los rendimientos
adf.test(rendimientos, alternative = "stationary")## Warning in adf.test(rendimientos, alternative = "stationary"): p-value smaller
## than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: rendimientos
## Dickey-Fuller = -11.99, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationaryEl p-valor es bajo (0.01), lo que indica que podemos rechazar la hipótesis nula.
El modelo ARIMA (Autoregressive Integrated Moving Average) es ampliamente utilizado para modelar series de tiempo estacionarias. Analizamos las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF) para identificar los órdenes del modelo.
# Graficar ACF y PACF
ggAcf(rendimientos, lag.max = 20) + ggtitle("Función de Autocorrelación (ACF)")
ggPacf(rendimientos, lag.max = 20) + ggtitle("Función de Autocorrelación Parcial (PACF)")
Las gráficas nos ayudan a identificar posibles valores de p y q para
el modelo ARIMA(p, d, q). En los rendimientos financieros, a menudo se
observa poca autocorrelación significativa. También podríamos hacer que
el algoritmo seleccione los valores para el modelo. Para esto utilizamos
la función auto.arima(), y seleccionará estos usando
criterios de información (AICc).
# Selección automática del modelo ARIMA
modelo_arima <- auto.arima(rendimientos, seasonal = FALSE)
# Resumen del modelo
summary(modelo_arima)## Series: rendimientos
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## -0.3563 0.2687 0.0011
## s.e. 0.1880 0.1926 0.0004
##
## sigma^2 = 0.000388: log likelihood = 4270.74
## AIC=-8533.49 AICc=-8533.46 BIC=-8511.73
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 7.872438e-07 0.01967924 0.01368463 NaN Inf 0.6399816 0.00343316El modelo seleccionado es un ARIMA(1,0,0), es decir, un modelo ARMA con un término autorregresivo pero sin media móvil, realmente, un modelo AR.
Verificamos si el modelo ajustado cumple con los supuestos necesarios.
# Graficar los residuos
checkresiduals(modelo_arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 404.26, df = 338, p-value = 0.007673
##
## Model df: 2. Total lags used: 340- Residuos Estandarizados: Deben comportarse como ruido blanco (✅).
- ACF de Residuos: No debe mostrar autocorrelación significativa (✅)..
- Prueba de Ljung-Box: Un p-valor alto indica que no hay autocorrelación en los residuos (💅).
Podríamos verificar la heteroscedasticidad. En series financieras, es común que exista esta estructura en los datos.
# Prueba de Arch
library(FinTS)##
## Attaching package: 'FinTS'## The following object is masked from 'package:forecast':
##
## Acf## The following object is masked from 'package:gt':
##
## sp500ArchTest(resid(modelo_arima))##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: resid(modelo_arima)
## Chi-squared = 252.28, df = 12, p-value < 2.2e-16La prueba es significativa,
Para capturar la heterocedasticidad, utilizamos modelos GARCH
(Generalized Autoregressive Conditional Heteroskedasticity). Utilizamos
el paquete rugarch para ajustar un modelo GARCH.
# Instalar y cargar el paquete rugarch
#install.packages("rugarch")
library(rugarch)## Loading required package: parallel##
## Attaching package: 'rugarch'## The following object is masked from 'package:purrr':
##
## reduce# Especificar el modelo ARIMA(1,0,0)-GARCH(1,1)
especificacion <- ugarchspec(mean.model = list(armaOrder = c(1,0)),
variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
distribution.model = "norm")
# Ajustar el modelo
modelo_garch <- ugarchfit(spec = especificacion, data = rendimientos)
# Resumen del modelo
show(modelo_garch)##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001757 0.000419 4.19379 0.000027
## ar1 -0.006134 0.026550 -0.23105 0.817278
## omega 0.000016 0.000008 2.06039 0.039362
## alpha1 0.103685 0.024233 4.27876 0.000019
## beta1 0.850074 0.023870 35.61232 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001757 0.000927 1.89578 0.057989
## ar1 -0.006134 0.027399 -0.22389 0.822841
## omega 0.000016 0.000042 0.39018 0.696405
## alpha1 0.103685 0.119539 0.86738 0.385737
## beta1 0.850074 0.099974 8.50292 0.000000
##
## LogLikelihood : 4449.133
##
## Information Criteria
## ------------------------------------
##
## Akaike -5.2222
## Bayes -5.2063
## Shibata -5.2223
## Hannan-Quinn -5.2163
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.004904 0.9442
## Lag[2*(p+q)+(p+q)-1][2] 0.078519 1.0000
## Lag[4*(p+q)+(p+q)-1][5] 0.916086 0.9586
## d.o.f=1
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.2805 0.5964
## Lag[2*(p+q)+(p+q)-1][5] 1.7351 0.6822
## Lag[4*(p+q)+(p+q)-1][9] 3.3195 0.7053
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.1956 0.500 2.000 0.6583
## ARCH Lag[5] 2.7576 1.440 1.667 0.3269
## ARCH Lag[7] 3.7495 2.315 1.543 0.3841
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 6.5406
## Individual Statistics:
## mu 0.50049
## ar1 0.41967
## omega 0.77602
## alpha1 0.12120
## beta1 0.06688
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.2167 0.8284
## Negative Sign Bias 1.2057 0.2281
## Positive Sign Bias 0.2567 0.7974
## Joint Effect 3.1504 0.3690
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 63.95 9.048e-07
## 2 30 65.11 1.372e-04
## 3 40 98.52 4.701e-07
## 4 50 94.00 1.166e-04
##
##
## Elapsed time : 0.06619382- Mu (\(μ\)): El término de la media es significativo (\(p < 0.01\)), lo que indica que hay un rendimiento promedio distinto de cero.
- AR(1): El coeficiente AR(1) no es significativo (\(p ≈ 0.42\)), lo que sugiere que el término autorregresivo puede no ser necesario.
- Omega (\(ω\)): Parámetro significativo que representa la varianza incondicional.
- Alpha1 (\(α₁\)) y Beta1 (\(β₁\)): Ambos son altamente significativos, confirmando la presencia de efectos ARCH y GARCH.
Las pruebas diagnósticas indican que no aparentan haber ni autocorrelación significativa en los residuos ni en los residuos al cuadrado, así como tampoco aparentan haber sesgos significativos en los residuos (signos). El modelo GARCH parece haber capturado adecuadamente la heteroscedasticidad. Los parámetros podrían ser algo inestables (verificar Nyblom)
# Graficar los residuos estandarizados
plot(modelo_garch, which = "all")##
## please wait...calculating quantiles...
Podría considerar un modelo más sencillo sin AR(1), y con colas más amplias en distribución t.
# Especificación de un modelo GARCH(1,1) sin término AR(1)
especificacion_simple <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
distribution.model = "norm")
modelo_garch_simple <- ugarchfit(spec = especificacion_simple, data = rendimientos)
# Modelo GARCH con distribución t de Student
especificacion_t <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
distribution.model = "std")
modelo_garch_t <- ugarchfit(spec = especificacion_t, data = rendimientos)
modelo_garch_t##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001646 0.000356 4.6267 4e-06
## omega 0.000014 0.000003 5.3728 0e+00
## alpha1 0.098932 0.006987 14.1604 0e+00
## beta1 0.867302 0.013832 62.7042 0e+00
## shape 4.858500 0.465941 10.4273 0e+00
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001646 0.000368 4.4719 0.000008
## omega 0.000014 0.000005 2.9430 0.003250
## alpha1 0.098932 0.018175 5.4432 0.000000
## beta1 0.867302 0.015897 54.5578 0.000000
## shape 4.858500 0.701722 6.9237 0.000000
##
## LogLikelihood : 4520.935
##
## Information Criteria
## ------------------------------------
##
## Akaike -5.3066
## Bayes -5.2906
## Shibata -5.3066
## Hannan-Quinn -5.3007
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.01439 0.9045
## Lag[2*(p+q)+(p+q)-1][2] 0.07550 0.9372
## Lag[4*(p+q)+(p+q)-1][5] 0.93053 0.8752
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.444 0.5052
## Lag[2*(p+q)+(p+q)-1][5] 1.891 0.6444
## Lag[4*(p+q)+(p+q)-1][9] 3.516 0.6716
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.1402 0.500 2.000 0.7081
## ARCH Lag[5] 2.8330 1.440 1.667 0.3149
## ARCH Lag[7] 3.8455 2.315 1.543 0.3694
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 15.2148
## Individual Statistics:
## mu 0.6105
## omega 3.5534
## alpha1 0.5540
## beta1 0.4539
## shape 0.8685
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.2684 0.7884
## Negative Sign Bias 1.2593 0.2081
## Positive Sign Bias 0.1642 0.8696
## Joint Effect 3.3647 0.3387
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 16.61 0.6160
## 2 30 28.45 0.4941
## 3 40 27.78 0.9100
## 4 50 32.84 0.9631
##
##
## Elapsed time : 0.0626092Finalmente, quiero visualizar usando el ARIMA(1,0,0) el rendimiento futuro de AAPL por 15 días.
Uso para esto la función forecast(): Esta función del
paquete forecast genera pronósticos basados en el modelo ARIMA ajustado
(modelo_arima). Con h=15 indico que quiero
saber los próximos 15 periodos (días). Graficamos ese pronóstico con
autoplot, y algo de personalización.
# Pronóstico de los próximos 15 días
pronostico_arima <- forecast(modelo_arima, h = 15)
# Graficar el pronóstico
autoplot(pronostico_arima) +
labs(title = "Pronóstico ARIMA de los Rendimientos de AAPL",
x = "Tiempo",
y = "Rendimiento") +
theme_minimal()
# Pronóstico de la volatilidad futura
pronostico_garch <- ugarchforecast(modelo_garch, n.ahead = 10)
# Extraer la volatilidad pronosticada
volatilidad_pronosticada <- sigma(pronostico_garch)
# Mostrar los resultados
volatilidad_pronosticada## 1974-08-30
## T+1 0.01410121
## T+2 0.01435463
## T+3 0.01459223
## T+4 0.01481529
## T+5 0.01502495
## T+6 0.01522223
## T+7 0.01540803
## T+8 0.01558318
## T+9 0.01574841
## T+10 0.01590440- Los pronósticos nos ayudan a entender las posibles tendencias futuras en los rendimientos y la volatilidad.
Conclusiones
- Análisis Exploratorio: Identificamos que los precios de las acciones no son estacionarios, pero los rendimientos logarítmicos sí lo son.
- Modelado ARIMA: Ajustamos un modelo ARIMA para capturar la dinámica en los rendimientos.
- Volatilidad y GARCH: Detectamos heterocedasticidad y utilizamos un modelo GARCH para modelar la volatilidad condicional.
- Pronósticos: Generamos pronósticos de rendimientos y volatilidad.
Es importante tener en cuenta que los mercados financieros son influenciados por muchos factores impredecibles, y ningún modelo puede capturar completamente esa complejidad. Por lo tanto, los pronósticos deben interpretarse con cautela y complementarse con análisis cualitativos y contextualizados.
Recomiendo aplicar estos métodos a diferentes activos financieros o índices para fortalecer la comprensión. Además de R, existen otras herramientas como Python y EViews que también son útiles para el análisis de series de tiempo. EViews está en declive según tengo entendido, pero sigue siendo utilizado en varios espacios.
10 Qué aprendimos en este taller
Aprendimos la existencia de varios modelos, maneras de representar gráficamente sus resultados, y evaluar su utilidad. He añadido elementos que aumentan la descripción de pedazos del taller, y que buscan corregir el orden en el que se trabajaron ejemplos.